AVISOの海面力学高度・海面地衡流データの作図

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

run.sh

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

#!/bin/sh

flist=$1

if [ ! -f $flist ]; then

echo error in $0: no such file, $flist

exit 1

fi

lonl=100.

lonr=180.

latb=0.

latt=60.

indir="../23.aviso/data/"

outdir1="output/h/"

outdir2="output/uv/"

para=para.txt

while read fname

do

in1=h/$fname

in2=uv/$(echo $fname | sed -e "s/h_/uv_/")

if [ ! -f ${indir}$in1 ]; then

echo error in $0: no such file, $in1

exit 1

fi

if [ ! -f ${indir}$in2 ]; then

echo error in $0: no such file, $in2

exit 1

fi

tmp1=$(basename $in1)

tmp2=$(echo $tmp1 | sed -e "s/dt_upd_global_merged_madt_//")

ofle=$(echo $tmp2 | sed -e "s/¥.nc/.txt/")

ofle2=$(echo $ofle | sed -e "s/h_/uv_/")

echo "&para" > $para

echo indir=¥"$indir¥" >> $para

echo infle=¥"$in1¥" >> $para

echo infl2=¥"$in2¥" >> $para

echo ofle=¥"${outdir1}$ofle¥" >> $para

echo ofle2=¥"${outdir2}$ofle2¥" >> $para

echo lonl=$lonl >> $para

echo lonr=$lonr >> $para

echo latb=$latb >> $para

echo latt=$latt >> $para

echo "&end" >> $para

# nc2txt < $para

plth.sh ${outdir1}$ofle

# plthuv.sh ${outdir1}$ofle ${outdir2}$ofle2

done < $flist

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

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=ifc

FC=gfortran

#FC=g95

#FC=g77

#FC=f90

#FC=f77

#FC=fort

#j ターゲット名(最終的に作りたい実行ファイル名)

TARGET1=nc2txt

TARGET2=

TARGETS=$(TARGET1) $(TARGET2)

#j OBJDIR (オブジェクトファイルをおくディレクトリ)

OBJDIR=obj

#j MOD (Fortran90以降のモジュールファイル)

MOD=

#j SRC7 (Fortran77のソースファイル)

SRC7=

#j SRC9 (Fortran90のソースファイル)

SRC9=nc2txt.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 コンパイルオプション(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 = #-convert big_endian -module obj

#-----------------------------------------

#j リンク用のオプション

#-----------------------------------------

LDFLAGS= -L/usr/local/lib -lnetcdf #-module ${OBJDIR}

#-----------------------------------------

#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

#j tarファイルを作る

tar:

tar cvf $(TARGET1).tar ./*.f90 ./*.f ./*.c ./*.h ./*.txt makefile

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

nc2txt.f90

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

program nc2txt

! Description:

!

! Author: aym

!

! Host: oceani18.fish.nagasaki-u.ac.jp

! Directory: /home/aym/11.Work10/96.AvisoMSLA-MADT/01.sample_data

!

! Revision history:

! 2010-10-17 16:03

! Initial Version

!

! References

! Tips: NetCDF. http://iprc.soest.hawaii.edu/users/sasakiyo/Japanese/misc/tips_netcdf.html

! NetCDF Tips. http://www.sci.hokudai.ac.jp/grp/poc/top/software/other/netcdf_tips/index.htm

! implicit none

! use

include 'netcdf.inc'

integer, parameter :: ysz=915 !721 ! 配列サイズ（緯度方向）

integer, parameter :: xsz=1080 !1440 ! 配列サイズ（経度方向）

integer :: ncid, status, rlatid, rlonid, rhid

character :: indir*100,infle*500,infl2*500,ofle*500,ofle2*500 !入出力ファイル名

character varname*10 !変数名

real,allocatable :: sla(:,:),u(:,:),v(:,:) ! データを収納するための変数

character :: latname*20,lonname*20

double precision,allocatable :: Lat(:), Lon(:) !緯度経度のデータ **倍精度**

real fillvalue !欠損値

real lonl,lonr,latb,latt

namelist /para/indir,infle,infl2,ofle,ofle2,lonl,lonr,latb,latt

write(*,'(a)')'Prog.: nc2txt'

read(5,para)

fillvalue=1.8E+19

! 配列寸法の決定

allocate(sla(ysz,xsz))

allocate(Lat(ysz), Lon(xsz))

varname='Grid_0001' ! slaの変数設定名

latname='NbLatitudes' !緯度

lonname='NbLongitudes' !経度

write(*,'(A,A,A)')' Input: ',trim(adjustl(indir)),trim(adjustl(infle))

write(*,*)

write(*,*)'xsz= ',xsz

write(*,*)'ysz= ',ysz

write(*,'(A,e12.4)')' Fillvalue (missing data) = ',fillvalue

write(*,*)

status = nf_open(trim(adjustl(indir))//trim(adjustl(infle)), nf_nowrite, ncid) ! ファイルのopenとNetCDF ID(ncid)の取得

if(status .ne. nf_noerr) call handle_err(status) !エラーチェック

! 緯度読み込み

write(*,'(a,a)')' Read ',trim(adjustl(latname))

status = nf_inq_varid(ncid, latname, rlatid) ! 変数のID取得

if(status .ne. nf_noerr) call handle_err(status)

status = nf_get_var_double(ncid, rlatid, lat)

if(status .ne. nf_noerr) call handle_err(status)

! 経度読み込み

write(*,'(a,a)')' Read ',trim(adjustl(lonname))

status = nf_inq_varid(ncid, lonname, rlonid) ! 変数のID取得

if(status .ne. nf_noerr) call handle_err(status)

status = nf_get_var_double(ncid, rlonid, lon)

if(status .ne. nf_noerr) call handle_err(status)

! sla読込み

write(*,'(a,a)')' Read ',trim(adjustl(varname))

status = nf_inq_varid(ncid, varname, rhid) ! 変数のID取得

if(status .ne. nf_noerr) call handle_err(status)

status = nf_get_var_real(ncid, rhid, sla)! 変数読み込み

if(status .ne. nf_noerr) call handle_err(status)

status = nf_close(ncid) ! ファイルのclose

if(status .ne. nf_noerr) call handle_err(status)

write(*,*)'Input file closed.'

write(*,*)''

! 配列寸法の決定

allocate(u(ysz,xsz),v(ysz,xsz))

write(*,'(A,A,A)')' Input: ',trim(adjustl(indir)),trim(adjustl(infl2))

status = nf_open(trim(adjustl(indir))//trim(adjustl(infl2)), nf_nowrite, ncid) ! ファイルのopenとNetCDF ID(ncid)の取得

if(status .ne. nf_noerr) call handle_err(status) !エラーチェック

! 緯度読み込み

write(*,'(a,a)')' Read ',trim(adjustl(latname))

status = nf_inq_varid(ncid, latname, rlatid) ! 変数のID取得

if(status .ne. nf_noerr) call handle_err(status)

status = nf_get_var_double(ncid, rlatid, lat)

if(status .ne. nf_noerr) call handle_err(status)

! 経度読み込み

write(*,'(a,a)')' Read ',trim(adjustl(lonname))

status = nf_inq_varid(ncid, lonname, rlonid) ! 変数のID取得

if(status .ne. nf_noerr) call handle_err(status)

status = nf_get_var_double(ncid, rlonid, lon)

if(status .ne. nf_noerr) call handle_err(status)

! u読込み

varname='Grid_0001' ! uの変数設定名

write(*,'(a,a)')' Read ',trim(adjustl(varname))

status = nf_inq_varid(ncid, varname, rhid) ! 変数のID取得

if(status .ne. nf_noerr) call handle_err(status)

status = nf_get_var_real(ncid, rhid, u)! 変数読み込み

if(status .ne. nf_noerr) call handle_err(status)

! v読込み

varname='Grid_0002' ! vの変数設定名

write(*,'(a,a)')' Read ',trim(adjustl(varname))

status = nf_inq_varid(ncid, varname, rhid) ! 変数のID取得

if(status .ne. nf_noerr) call handle_err(status)

status = nf_get_var_real(ncid, rhid, v)! 変数読み込み

if(status .ne. nf_noerr) call handle_err(status)

status = nf_close(ncid) ! ファイルのclose

if(status .ne. nf_noerr) call handle_err(status)

write(*,*)'Input file closed.'

write(*,*)''

! 出力

call system( &

& 'if [ ! -d output/h ]; then &

& (echo "Create output directory"; mkdir -p output/h ); &

& fi')

call system( &

& 'if [ ! -d output/uv ]; then &

& (echo "Create output directory"; mkdir -p output/uv ); &

& fi')

iu1=21

write(*,'(A,A)')' Output: ',trim(adjustl(ofle))

open(iu1,file=ofle)! ,form="unformatted")

write(*,'(A,A,A)')' Open ',trim(adjustl(ofle)),'.'

! call print_header(iu2)

iu2=22

write(*,'(A,A)')' Output: ',trim(adjustl(ofle2))

open(iu2,file=ofle2)! ,form="unformatted")

write(*,'(A,A,A)')' Open ',trim(adjustl(ofle2)),'.'

! slaの1番目の添え字は緯度を表し、2番目の添え字は経度を表している

do j=1,ysz

do i=1,xsz

if(sla(j,i).le.fillvalue)then

if(lon(i)>=lonl.and.lon(i)<=lonr.and.lat(j)>=latb.and.lat(j)<=latt)then

write(iu1,'(f11.5,f10.5,f10.2)')lon(i),lat(j),sla(j,i)

endif

endif

enddo !i

enddo !j

close(iu1)

do j=1,ysz

do i=1,xsz

if(u(j,i).le.fillvalue)then

if(lon(i)>=lonl.and.lon(i)<=lonr.and.lat(j)>=latb.and.lat(j)<=latt)then

write(iu2,'(f11.5,f10.5,2f11.3)')lon(i),lat(j),u(j,i),v(j,i)

endif

endif

enddo !i

enddo !j

close(iu2)

write(*,'(a)')'Done.'

write(*,*)

stop

end program nc2txt

subroutine handle_err(status)

integer,intent(in) :: status

write(*,*)nf_strerror(status)

stop

end subroutine handle_err

subroutine print_header(iunit)

! Description: ヘッダ情報の書き出し

! このサブルーチンは, gfortran, ifortで動作することを確認している.

!

! 作成したデータが、何時、何処で、誰が、どのプログラムを使って作成した

! か記録しておくと、後々作業をし直す必要が生じたときに役に立つ.

! また、以前使ったプログラムを一部改変して利用するときなど、そのプログ

! ラムを探すのに役に立つことが多い

!

!

! Author: aym

!

integer,intent(in) :: iunit

!

character dirinfo*144, userinfo*72, hostinfo*72

INTEGER date_time(8)

CHARACTER REAL_CLOCK(3)*12

character progname*100

integer idxpn, idx

character shellinfo*80

!

! Get date, time, current directory, username, and hostname

! p.345 of for_lang.pdf (intel fortran user manual)

! DATE_AND_TIME is an intrisic S/R and works on most of Fotran 90

! Compilers. g77 also supports this S/R.

CALL DATE_AND_TIME (REAL_CLOCK (1), REAL_CLOCK (2), &

& REAL_CLOCK (3), date_time)

call getenv("PWD", dirinfo)

call getenv("USER", userinfo)

!j 使用しているシェルのチェック

!j 現在, sh, bash, csh, tcshのみに対応

call getenv("SHELL",shellinfo)

idxcsh=index(shellinfo,"/csh")

idxtcsh=index(shellinfo,"/tcsh")

idxsh=index(shellinfo,"/sh")

idxbash=index(shellinfo,"/bash")

if(idxcsh.ne.0.or.idxtcsh.ne.0)then

call getenv("HOST", hostinfo) !csh, tcsh, etc.

else if(idxsh.ne.0.or.idxbash.ne.0)then

call getenv("HOSTNAME", hostinfo) !sh, bash, etc.

endif

write(iunit,'(a, &

& i2.2,a,i2.2,a,i4.4, a,i2.2,a,i2.2,a,i2.2, a,i3.2,a)') &

& '# Date and time: ', &

& date_time(2),'/',date_time(3),'/',date_time(1), &

& ' at ',date_time(5),':',date_time(6),':',date_time(7), &

& ' ',-date_time(4)/60,':00'

write(iunit,'(A,A)')'# hostname: ',trim(adjustl(hostinfo))

write(iunit,'(A,A)')'# cwd: ',trim(adjustl(dirinfo))

write(iunit,'(A,A)')'# user: ',trim(adjustl(userinfo))

! Get program name

call getarg( 0, progname)

write(iunit,'(A,A)')'# program name: ',trim(adjustl(progname))

return

end subroutine print_header

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

para.txt

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

&para

indir="../23.aviso/data/"

infle="h/dt_upd_global_merged_madt_h_20091230_20091230_20100921.nc"

infl2="uv/dt_upd_global_merged_madt_uv_20091230_20091230_20100921.nc"

ofle="output/h/h_20091230_20091230_20100921.txt"

ofle2="output/uv/uv_20091230_20091230_20100921.txt"

lonl=100.

lonr=180.

latb=0.

latt=60.

&end

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

plthuv.sh

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

#!/bin/sh

rm -rf .gmtcommand* .gmtdefaults*

#

# 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="¥312" # degree symbol

# 色のRGB値を設定 (線や点に色をつけるときRGB値を直接指定するより分かりやすい)

red="255/0/0"

green="0/255/0"

blue="0/0/255"

# 数字の出力フォーマット(project, grd2xzy等で使う)

gmtset D_FORMAT %lg #デフォルト(規定値)

# 桁数を指定（例：123.45678)

#gmtset D_FORMAT %12.6f

# 状況に応じて

#gmtset D_FORMAT %12.5g

#

# 引数のチェック

#

if [ $# == 0 ]; then

echo Error in $0 : No argument

echo Usage: $0 FILENAME

echo

exit 1

fi

# 最初の引数を入力ファイル名とする．

in1=$1

in2=$2

if [ ! -f $in1 ]; then

echo Error in $0 : No input file, $in1

exit 1

fi

if [ ! -f $in2 ]; then

echo Error in $0 : No input file, $in2

exit 1

fi

# 入力ファイルの拡張子を指定

# 拡張子が3文字であることを仮定しているので注意

ext="txt"

# 入力ファイル名から拡張子を取得

# 拡張子が3文字であることを仮定しているので注意

ext_infile=${in1:${#in1}-3:${#in1}}

# 拡張子のチェック

if [ ${ext_infile} != ${ext} ]; then

echo "Error in $0 : Extension of the input file name must be ${ext}."

exit 1

fi

outdir="ps/huv/"

if [ ! -d $outdir ]; then

mkdir -p $outdir

fi

tmp=$(echo $in2 | sed -e "s/uv_/huv_/")

out=${outdir}$(basename $tmp .${ext}).ps

# 入力ファイルを上書きしないよう，念のためチェック

if [ $in1 == $out ]; then

echo Error in $0 : Input file name is identical to output file name.

exit 1

fi

echo Input:

echo $in1

echo $in2

echo Output: $out

# ---------------------------------------------------------------------------------

# ここから下にgmtのコマンドを書く

# ---------------------------------------------------------------------------------

grd=$0.tmp.$$.grd

grd2=$0.tmp2.$$.grd

grd3=$0.tmp3.$$.grd

cpt=$0.$$.cpt

reso=0.2/0.2 # Resolution of the plot

range=120/160/15/48 # Plot range

maptype=M

size=5

xanot=a10f5

yanot=a10f5

makecpt -Cjet -T-100/100/5 > $cpt

echo

echo "Now working. Be patient..."

echo

awk '{if ($3 !="NaN" && $1 !="#") printf "%13.7f %12.7f %9.2f¥n", $1, $2, $3}' ${in1} |¥

surface -R${range} -T1 -I${reso} -G${grd}

grdgradient ${grd} -G${grd2} -D -S${grd3}

grdimage ${grd3} -C$cpt -R${range} -J${maptype}${size} -P -K -X1.5 -Y2 > $out

grdcontour ${grd} -R -J${maptype} -A30f10 -C10 -L30/1000 -W1 ¥

-O -K >> $out

grdcontour ${grd} -R -J${maptype} -A30f10 -C10 -L-1000/-30 -W1t15_5:15 -O -K >> $out

# Vector Size

VSIZE=100 #0.1 # cms-1

VSCALE=0.1 # inch

VFACTOR=`echo "scale=6; ${VSCALE}/${VSIZE}" | bc`

# Vectors with color

#awk -v SC="${VFACTOR}" '{if ( ($1 !="#") && ((NR % 4)==0) ) ¥

# printf "%13.7f %12.7f %10.3f %10.3f %10.4f¥n", ¥

# $1, $2, sqrt($3*$3+$4*$4),(180.0/3.1415926535)*atan2($4, $3), SC*sqrt($3*$3+$4*$4)}' ${in} |¥

# psxy -R${range} -J${maptype}${size} -Sv0.01/0.1/0.03n0.2 -C$cpt ¥

# -O -K >> $out

# Black and While

awk -v SC="${VFACTOR}" '{if ( ($1 !="#") && ((NR % 2)==0) ) ¥

printf "%13.7f %12.7f %10.3f %10.4f¥n", ¥

$1, $2, (180.0/3.1415926535)*atan2($4, $3), SC*sqrt($3*$3+$4*$4)}' ${in2} |¥

psxy -R${range} -J${maptype}${size} -Sv0.01/0.1/0.03n0.2 -G0 ¥

-O -K >> $out

pscoast -R$range -J${maptype}${size} -G180 -W3 -Dl -B${xanot}/${yanot}WSne ¥

-L125/43/43/300 -O -K >> $out

echo "pscoast command finished."

tmp=$(basename $in1)

day="${tmp:2:4}/${tmp:6:2}/${tmp:8:2}"

pstext << end -R$range -J${maptype}${size} -O -K >> $out

122 46 18 0 1 LM $day

end

psscale -D5.5/2/4/0.1 -C$cpt -Ba20f10:""::."cm/s": -O -K ¥

>> $out

rm -rf $grd $grd2 $grd3 $cpt

echo "Delete temp. files"

# ---------------------------------------------------------------------------------

# gmtのコマンドはここまで

# ---------------------------------------------------------------------------------

# 図を作成した時刻，現在のディレクトリ名，入出力ファイル名，コメントなどを図に入れる．

comment="$(ls -l $in1 | awk '{print $6, $7, $8}')" # コメント

# yoffset: y方向（紙の縦方向）のオフセット．

# 単位はインチ（1インチ=2.54cm, A4の縦はおおよそ11インチ）．

yoffset=5.3

currentdir=$(pwd)

if [ ${#currentdir} -gt 90 ]; then

curdir1=${currentdir:1:90}

curdir2=${currentdir:91}

else

curdir1=${currentdir}

curdir2="¥ "

fi

now=$(date)

host=$(hostname)

pstext -JX6/1.2 -R0/1/0/1.2 -N -Y${yoffset} << 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} (${comment})

0 0.1 9 0 1 LM Output: ${out}

EOF

echo "$0 finished."

echo

exit 0