ジオポテンシャル高度

------------------------------

List of the following files:

------------------------------

run.sh

pl_cont_pstereo.sh

gmtpar.sh

note.sh

src/caldat.f90

src/julday.f90

src/read_ncep_4d.f90

src/read_ncep_4d_int.f90

src/makefile

------------------------------

Machine info

------------------------------

localhost

/work2/am/TEACHING/Global.Meteorology/GPH

Wed Nov 1 11:04:47 JST 2017

======================

run.sh

======================

#!/bin/bash

# Description:

#

# Author: am

#

# Host: aofd30

# Directory: /work2/am/teaching/JetStream_Contours/plot

#

# Revision history:

# This file is created by /usr/local/mybin/nbscr.sh at 15:21 on 11-04-2011.

echo "Shell script, $(basename $0) starts."

echo

if [ $# -ne 3 ]; then

echo Error in $0 : Wrong number of arguments.

echo "Usage : $0 <varname> <yyyymmdd1> <yyyymmdd1> <inthr>"

echo "varname: variable name (e.g., hgt)"

echo "yyyymmdd1 : year month day (e.g., 20120920)"

echo "yyyymmdd2 : year month day (e.g., 20121020)"

exit 1

fi

varname=$1

yyyymmdd1=$2

yyyymmdd2=$3

yyyy1=${yyyymmdd1:0:4}

mm1=${yyyymmdd1:4:2}

dd1=${yyyymmdd1:6:2}

yyyy2=${yyyymmdd2:0:4}

mm2=${yyyymmdd2:4:2}

dd2=${yyyymmdd2:6:2}

if [ $yyyy1 -ne $yyyy2 ]; then

echo Error in $0 : year must be the same

echo yyyymmdd1 $yyyymmdd1

echo yyyymmdd2 $yyyymmdd2

exit 1

fi

outdir="output/"

if [ ! -d $outdir ]; then

mkdir -p $outdir

fi

indir=/misc/raid112MANDA/data/NCEP1

k=8 # Vertical level

im=144 # See header of the input file (netCDF file). Type "ncdump -h INPUT_FILE_NAME"

jm=73

km=17

nm=1204

#nm=1204

exe=./read_ncep_4d

if [ ! -f $exe ]; then

echo Error in $0 : No such file, $exe

exit 1

fi

if [ ! -f $infle ]; then

echo Error in $0 : No such file, $infle

exit 1

fi

infle="${indir}/${varname}.${yyyy1}.nc"

echo

echo $infle

echo

jday1=$(date -d${yyyy1}${mm1}${dd1} +%j)

jday2=$(date -d${yyyy2}${mm2}${dd2} +%j)

nday=$(expr $jday2 - $jday1 )

jday0=$(date -d"${yyyy1}/1/1" '+%j')

export LANG=C

iday=0

while [ $iday -le ${nday} ]; do

yyyymmdd=$(date -d"${yyyy1}/${mm1}/${dd1} ${iday} day" '+%Y%m%d')

jday=$(date -d"${yyyymmdd}" '+%j')

jday=$(echo "scale=0; $jday - 1"|bc)

imh=4

i=1

while [ $i -le $imh ]; do

hh=$(echo "scale=0; ($i - 1) * 6" |bc)

nstep=$(echo "scale=0; (${jday} - ${jday0})*4 + $i"|bc)

echo "yyyymmdd=${yyyymmdd} hh=${hh}"

date -d"${yyyy1}/1/1 ${jday} day" '+%Y%m%d'

echo "jday0=${jday0} jday=${jday} nstep=${nstep}"

echo

cat << EOF > para.txt

&para

varname="${varname}",

yyyymmdd="${yyyymmdd}",

hh=${hh},

infle="${infle}",

nstep=${nstep},

klevel=${k},

im=${im},

jm=${jm},

km=${km},

nm=${nm},

&end

EOF

$exe < para.txt

i=$(expr $i + 1)

done #i

iday=$(expr $iday + 1)

done #jday

echo "Done $(basename $0)"

echo

----------------------

End of run.sh

----------------------

======================

pl_cont_pstereo.sh

======================

#!/bin/bash

# Description:

#

# Author: am

#

# Host: aofd30

# Directory: /work2/am/teaching/JetStream_Contours/plot

#

# Revision history:

# This file is created by /usr/local/mybin/ngmt.sh at 13:41 on 11-04-2011.

. ./gmtpar.sh

export LANG=C

echo "Bash script $0 starts."

range=0/360/20/85

size=0/90/1.27i/60

xanot=a30g10

yanot=a30g10

anot=${xanot}/${yanot}WSe

#indir="../output/"

if [ $# -ne 1 ]; then

echo "Error in $0 : Wrong number of arguments."

echo Usage: $0 input #number

echo input : input file name

# echo number : file number

exit

fi

in=$1 #${indir}$1

#no=$2

if [ ! -f $in ]; then

echo Error in $0 : No such file, $in

exit 1

fi

psdir="./Fig/cont/"

if [ ! -d ${psdir} ];then

mkdir -p $psdir

fi

out=${psdir}$(basename $in .asc)_cont.ps

echo Input : $in

echo Output : $out

xoffset=1.5

yoffset=12.5

comment=

. ./note.sh

grd=$0.grd

reso=2.5/2.5

awk '{if ($1 != "#") print $1,$2,$3/10}' $in | \

surface -R$range -I${reso:-"10/5"} -T1 -G$grd

yyyymmdd=$(awk '{if ($2 == "yyyymmdd,hr=") print $3}' $in)

hr=$(awk '{if ($2 == "yyyymmdd,hr=") print $4}' $in)

reftime=$(awk '{if ($2 == "yyyymmdd,hr=") print $5}' $in)

title="${hr}${reftime}${yyyymmdd}"

pscoast -R$range -Js$size \

-B$anot:.${title}: -Dc -G255 -W1 \

-X-.5 -Y-9.5 -O -K >> $out

grdcontour $grd -R$range -Js$size -C10 -A20f10 -W2/150 \

-O -K \

>> $out

range=0/1/0/1

size=X5/11

pstext << end -R$range -J$size -O -K -P -Y4 >> $out

0.1 0.4 20 0 0 LM Gakka:

0.6 0.4 20 0 0 LM Program:

0.1 0.35 20 0 0 LM Bango:

0.6 0.35 20 0 0 LM Namae:

end

rm $grd

echo "Done $0.sh"

----------------------

End of pl_cont_pstereo.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"

yellow="255/255/0"

purple="160/32/240"

magenta="255/0/255"

white="255/255/255"

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

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

# 桁数を指定（例：123.45678)

#gmtset D_FORMAT %12.6f

# 状況に応じて

#gmtset D_FORMAT %12.5g

----------------------

End of gmtpar.sh

----------------------

======================

note.sh

======================

# Print time, current working directory and output filename

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 \

-K -P > $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}

----------------------

End of note.sh

----------------------

======================

src/caldat.f90

======================

subroutine caldat(JULIAN,IYYY,MONTH,DD)

! Description:

!

! Author: am

!

! Host: aofd30

! Directory: /work2/am/12.Work11/75.Mitsui/52.Vertical_Profile/src

!

! Revision history:

! This file is created by /usr/local/mybin/nff.sh at 11:00 on 11-01-2011.

! use

! implicit none

! write(*,'(a)')'Subroutine: caldat'

! write(*,*)''

!***********************************************************************

!

! ********** SUBROUTINE DESCRIPTION:

!

! GIVEN THE JULIAN DAY, RETURNS THE YEAR, MONTH AND DAY OF MONTH.

!

! ********** ORIGINAL AUTHOR AND DATE:

!

! PRESS,FLANNERY,TEUKOLSKY,VETTERLING 1986.

! NUMERICAL RECIPES

!

! ********** REVISION HISTORY:

!

!

! ********** ARGUMENT DEFINITIONS:

!

INTEGER,intent(in):: JULIAN

integer,intent(out)::IYYY,MONTH,DD

!

! NAME IN/OUT DESCRIPTION

!

! JULIAN I THE JULIAN DAY

! IYYY O THE YEAR

! MONTH O THE MONTH (1 TO 12)

! DD O THE DAY OF THE MONTH

!

! ********** COMMON BLOCKS:

!

! NONE

!

! ********** LOCAL PARAMETER DEFINITIONS:

!

INTEGER IGREG

PARAMETER (IGREG=2299161)

!

! ********** LOCAL VARIABLE DEFINITIONS:

!

INTEGER JALPHA,JA,JB,JC,JD,JE

!

! NAME DESCRIPTION

!

!

! ********** OTHER ROUTINES AND FUNCTIONS CALLED:

!

!

!---+67--1----+----2----+----3----+----4----+----5----+----6----+----7--

!

!

IF (JULIAN .GE. IGREG) THEN

JALPHA = INT(((JULIAN - 1867216) - 0.25)/36524.25)

JA = JULIAN + 1 + JALPHA - INT(0.25*JALPHA)

ELSE

JA = JULIAN

ENDIF

JB = JA + 1524

JC = INT(6680. + ((JB - 2439870) - 122.1)/365.25)

JD = 365*JC + INT(0.25*JC)

JE = INT((JB - JD)/30.6001)

DD = JB - JD - INT(30.6001*JE)

MONTH = JE - 1

IF (MONTH .GT. 12) MONTH = MONTH - 12

IYYY = JC - 4715

IF (MONTH .GT. 2) IYYY = IYYY - 1

IF (IYYY .LE. 0) IYYY = IYYY - 1

RETURN

! write(*,'(a)')'Done subroutine caldat.'

! write(*,*)

end subroutine caldat

----------------------

End of src/caldat.f90

----------------------

======================

src/julday.f90

======================

INTEGER FUNCTION JULDAY(IYYY,MONTH,DD)

! Description:

!

! Author: am

!

! Host: aofd30

! Directory: /work2/am/12.Work11/75.Mitsui/52.Vertical_Profile/src

!

! Revision history:

! This file is created by /usr/local/mybin/nff.sh at 10:59 on 11-01-2011.

!

! ********** SUBROUTINE DESCRIPTION:

!

! FINDS THE JULIAN DAY FROM A DATE.

!

! ********** ORIGINAL AUTHOR AND DATE:

!

! PRESS,FLANNERY,TEUKOLSKY,VETTERLING 1986.

! NUMERICAL RECIPES

!

! ********** REVISION HISTORY:

!

!

! ********** ARGUMENT DEFINITIONS:

!

INTEGER,intent(inout):: IYYY,MONTH,DD

!

! NAME IN/OUT DESCRIPTION

!

! IYYY I YEAR

! MONTH I MONTH (1 TO 12)

! DD I DAY OF MONTH

! JULDAY O JULIAN DAY

!

! ********** COMMON BLOCKS:

!

! NONE

!

! ********** LOCAL PARAMETER DEFINITIONS:

!

INTEGER IGREG

PARAMETER (IGREG = 15 + 31*(10 + 12*1582))

!

! ********** LOCAL VARIABLE DEFINITIONS:

!

INTEGER JY,JM,JA

!

! NAME DESCRIPTION

!

!

! ********** OTHER ROUTINES AND FUNCTIONS CALLED:

!

! INT - INTRINSIC TRUNCATE

!

!---+67--1----+----2----+----3----+----4----+----5----+----6----+----7--

!

IF (IYYY .LT. 0) IYYY = IYYY + 1

IF (MONTH .GT. 2) THEN

JY = IYYY

JM = MONTH + 1

ELSE

JY = IYYY - 1

JM = MONTH + 13

ENDIF

JULDAY = INT(365.25*JY) + INT(30.6001*JM) + DD + 1720995

IF (DD + 31*(MONTH + 12*IYYY) .GE. IGREG) THEN

JA = INT(0.01*JY)

JULDAY = JULDAY + 2 - JA + INT(0.25*JA)

ENDIF

RETURN

end function julday

----------------------

End of src/julday.f90

----------------------

======================

src/read_ncep_4d.f90

======================

program read_ncep_4d

! Description:

!

! Author: am

!

! Host: aofd30

! Directory: /work2/am/teaching/JetStream_Contours/src

!

! Revision history:

! This file is created by /usr/local/mybin/nff.sh at 13:23 on 11-04-2011.

! use

! implicit none

! include '/usr/local/include/netcdf.inc'

! include '/usr/local/netcdf-4.1.3/include/netcdf.h'

include '/usr/local/netcdf-4.1.3/include/netcdf.inc'

! ncid:ファイルのID番号、 varid: 変 数のID番号

integer ncid,varid

character varname*10

character infle*500,ofle*500

character var*50

real,allocatable :: level(:),lon(:),lat(:)

real(8),allocatable :: time(:)

! integer(2),allocatable :: sdata(:,:,:,:)

real,allocatable :: sdata(:,:,:,:)

integer,allocatable :: dimids(:),nshape(:)

character*70 err_message

character*100 dimname !dimnameは各次元の名前。

character(len=8):: yyyymmdd

integer hh

namelist /para/varname,yyyymmdd,hh,infle,nstep,klevel,im,jm,km,nm

write(*,'(a)')'Program read_ncep_4d starts.'

write(*,*)''

read(*,nml=para)

print *,'varname=',trim(varname)

print *,'yyyymmdd=',yyyymmdd

print *,'hh=',hh

print *,'infle=',trim(infle)

print *,'nstep=',nstep

print *,'klevel=',klevel

print *,'im=',im

print *,'jm=',jm

print *,'km=',km

print *,'nm=',nm

allocate(level(km),lon(im),lat(jm),time(nm),sdata(im,jm,km,nm))

print '(A,A)','Input: ',infle(1:lnblnk(infle))

sta = nf_open(infle, nf_nowrite, ncid) ! ファイルのopenとNetCDF ID(ncid)の取得

print *,'ncid= ',ncid

print *,'status= ',sta

print *

!もし、openが上手く行かなかったら、sta /= 0となる。

if(sta.ne.0) then

! その時のエラー を教えてくれる。

print *,'ERROR : Can NOT open file.'

! err_message = nf_strerror(sta)

! write(6,*) err_message

stop

endif

var=varname(1:lnblnk(varname))

print '(A)',var(1:lnblnk(var))

! ファイルID(ncid)からvarで設定した変数のID(varid)を得る。

sta = nf_inq_varid(ncid,var,varid)

print *,'varid=', varid

print *,'sta= ', sta

print *

! スケールファクターを得る。(*1)

! sta = nf_get_att_real(ncid,varid,'scale_factor',scale)

! print *,'scale= ',scale

! print *,'sta= ',sta

! print *

! オフセットの 値を得る(*2)

! sta = nf_get_att_real(ncid,varid,'add_offset',offset)

! print *,'offset=', offset

! print *,'sta= ',sta

! print *

! varidからndimsを得る。ndimsとは次元の数。二次元データなら2

stat = nf_inq_varndims(ncid, varid, ndims)

print *,'ndims= ', ndims

print *

! varidからdimidsを得る。dimidsとはそれぞれの次元に対するID

allocate(dimids(ndims))

stat = nf_inq_vardimid(ncid, varid, dimids)

do i=1,ndims

write(6,'("dimids(",i2,")=",i9 )') i,dimids(i)

enddo

print *

allocate(nshape(ndims))

do i=1,ndims

! dimidsからnshapeを得る。

! nshapeとはそれぞれの次元に対する最大データ数 (格子数)

stat = nf_inq_dimlen(ncid, dimids(i), nshape(i))

write(6,'("nshape(",i2,")=",i9 )') i,nshape(i)

enddo !i

print *

do i=1,ndims

stat = nf_inq_dimname(ncid, dimids(i), dimname)

write(6,'("dimname(",i2,")=",1x,a23)') i,dimname

enddo

print *

print *,'Read ',var(1:lnblnk(var))

! stat = nf_get_var_int2(ncid, varid, sdata)

stat = nf_get_var_real(ncid, varid, sdata)

print *,'stat= ',stat

!

var='lon'

print '(A)',var(1:lnblnk(var))

! ファイルID(ncid)からvarで設定した変数のID(varid)を得る。

sta = nf_inq_varid(ncid,var,varid)

print *,'varid=', varid

print *,'sta= ', sta

print *,'Read ',var(1:lnblnk(var))

stat = nf_get_var_real(ncid, varid, lon)

print *,'stat= ',stat

!

var='lat'

print '(A)',var(1:lnblnk(var))

! ファイルID(ncid)からvarで設定した変数のID(varid)を得る。

sta = nf_inq_varid(ncid,var,varid)

print *,'varid=', varid

print *,'sta= ', sta

print *,'Read ',var(1:lnblnk(var))

stat = nf_get_var_real(ncid, varid, lat)

print *,'stat= ',stat

print *

var='time'

print '(A)',var(1:lnblnk(var))

! ファイルID(ncid)からvarで設定した変数のID(varid)を得る。

sta = nf_inq_varid(ncid,var,varid)

print *,'varid=', varid

print *,'sta= ', sta

print *,'Read ',var(1:lnblnk(var))

stat = nf_get_var_double(ncid, varid, time)

print *,'stat= ',stat

print *

var='level'

print '(A)',var(1:lnblnk(var))

! ファイルID(ncid)からvarで設定した変数のID(varid)を得る。

sta = nf_inq_varid(ncid,var,varid)

print *,'varid=', varid

print *,'sta= ', sta

print *,'Read ',var(1:lnblnk(var))

stat = nf_get_var_real(ncid, varid, level)

print *,'stat= ',stat

print *

print *,'klevel=',klevel

print *,'level=',level(klevel),' [hPa]'

print *

! Note: time index (n) is fixed and set to 1 here.

print '(A,A8,1x,i2.2)','yyyymmdd,hh= ',yyyymmdd,hh

is=1;ie=len("output/")

write(ofle(is:ie),'(A)')"output/"

is=ie+1;ie=is+lnblnk(varname)-1

write(ofle(is:ie),'(A)')varname(1:lnblnk(varname))

ilevel=int(level(klevel))

is=ie+1;ie=is+3

write(ofle(is:ie),'(i3.3,A)')ilevel,'_'

is=ie+1;ie=is+10

write(ofle(is:ie),'(A8,A,i2.2)')yyyymmdd,'_',hh

is=ie+1;ie=is+len(".asc")-1

write(ofle(is:ie),'(A)')".asc"

print '(A)',ofle(1:lnblnk(ofle))

print '(A,A)','Output: ',ofle(1:lnblnk(ofle))

print *

open(20,file=ofle)

write(20,'(A,A)') '# input= ',infle(1:lnblnk(infle))

write(20,'(A,A8,A,i2.2,A)')'# yyyymmdd,hr= ',&

& yyyymmdd,' ',hh,' Z'

write(20,'(A,f10.2,A)')'# level= ',level(klevel),' hPa'

do j=1,jm

do i=1,im

! varout=sdata(i,j,klevel,nstep)*scale + offset

varout=sdata(i,j,klevel,nstep)

write(20,'(2f12.6,f12.4,1x,2i5)') &

lon(i),lat(j),varout,i,j

enddo !i

enddo !j

close(20)

! *1： scaleとは、ある変数の平均的なオーダー

! *2： offsetとは、ある変数の平均的な値

! data(i,j) = sdata(i,j) * scale + offset

! のようにして、本当のデータを得る

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

write(*,*)

end program read_ncep_4d

----------------------

End of src/read_ncep_4d.f90

----------------------

======================

src/read_ncep_4d_int.f90

======================

program read_ncep_4d

! Description:

!

! Author: am

!

! Host: aofd30

! Directory: /work2/am/teaching/JetStream_Contours/src

!

! Revision history:

! This file is created by /usr/local/mybin/nff.sh at 13:23 on 11-04-2011.

! use

! implicit none

include '/usr/local/include/netcdf.inc'

! ncid:ファイルのID番号、 varid: 変 数のID番号

integer ncid,varid

character varname*10

character infle*500,ofle*500

character var*50

real,allocatable :: level(:),lon(:),lat(:)

real(8),allocatable :: time(:)

integer(2),allocatable :: sdata(:,:,:,:)

integer,allocatable :: dimids(:),nshape(:)

character*70 err_message

character*100 dimname !dimnameは各次元の名前。

namelist /para/varname,infle,nstep,klevel,im,jm,km,nm

write(*,'(a)')'Program read_ncep_4d starts.'

write(*,*)''

open(1,file="para.txt",action="read")

read(1,para)

close(1)

allocate(level(km),lon(im),lat(jm),time(nm),sdata(im,jm,km,nm))

print '(A,A)','Input: ',infle(1:lnblnk(infle))

sta = nf_open(infle, nf_nowrite, ncid) ! ファイルのopenとNetCDF ID(ncid)の取得

print *,'ncid= ',ncid

print *,'status= ',sta

print *

!もし、openが上手く行かなかったら、sta /= 0となる。

if(sta.ne.0) then

! その時のエラー を教えてくれる。

print *,'ERROR : Can NOT open file.'

err_message = nf_strerror(sta)

write(6,*) err_message

stop

endif

var=varname(1:lnblnk(varname))

print '(A)',var(1:lnblnk(var))

! ファイルID(ncid)からvarで設定した変数のID(varid)を得る。

sta = nf_inq_varid(ncid,var,varid)

print *,'varid=', varid

print *,'sta= ', sta

print *

! スケールファクターを得る。(*1)

sta = nf_get_att_real(ncid,varid,'scale_factor',scale)

print *,'scale= ',scale

print *,'sta= ',sta

print *

! オフセットの 値を得る(*2)

sta = nf_get_att_real(ncid,varid,'add_offset',offset)

print *,'offset=', offset

print *,'sta= ',sta

print *

! varidからndimsを得る。ndimsとは次元の数。二次元データなら2

stat = nf_inq_varndims(ncid, varid, ndims)

print *,'ndims= ', ndims

print *

! varidからdimidsを得る。dimidsとはそれぞれの次元に対するID

allocate(dimids(ndims))

stat = nf_inq_vardimid(ncid, varid, dimids)

do i=1,ndims

write(6,'("dimids(",i2,")=",i9 )') i,dimids(i)

enddo

print *

allocate(nshape(ndims))

do i=1,ndims

! dimidsからnshapeを得る。

! nshapeとはそれぞれの次元に対する最大データ数 (格子数)

stat = nf_inq_dimlen(ncid, dimids(i), nshape(i))

write(6,'("nshape(",i2,")=",i9 )') i,nshape(i)

enddo !i

print *

do i=1,ndims

stat = nf_inq_dimname(ncid, dimids(i), dimname)

write(6,'("dimname(",i2,")=",1x,a23)') i,dimname

enddo

print *

print *,'Read ',var(1:lnblnk(var))

stat = nf_get_var_int2(ncid, varid, sdata)

print *,'stat= ',stat

!

var='lon'

print '(A)',var(1:lnblnk(var))

! ファイルID(ncid)からvarで設定した変数のID(varid)を得る。

sta = nf_inq_varid(ncid,var,varid)

print *,'varid=', varid

print *,'sta= ', sta

print *,'Read ',var(1:lnblnk(var))

stat = nf_get_var_real(ncid, varid, lon)

print *,'stat= ',stat

!

var='lat'

print '(A)',var(1:lnblnk(var))

! ファイルID(ncid)からvarで設定した変数のID(varid)を得る。

sta = nf_inq_varid(ncid,var,varid)

print *,'varid=', varid

print *,'sta= ', sta

print *,'Read ',var(1:lnblnk(var))

stat = nf_get_var_real(ncid, varid, lat)

print *,'stat= ',stat

print *

var='time'

print '(A)',var(1:lnblnk(var))

! ファイルID(ncid)からvarで設定した変数のID(varid)を得る。

sta = nf_inq_varid(ncid,var,varid)

print *,'varid=', varid

print *,'sta= ', sta

print *,'Read ',var(1:lnblnk(var))

stat = nf_get_var_double(ncid, varid, time)

print *,'stat= ',stat

print *

! Note: time index (n) is fixed and set to 1 here.

idate_0=julday(1,1,1)

date_n=time(nstep)/24.0 + idate_0

idate_n=int(date_n)

! print *,idate_0

! print *,idate_n

call caldat(idate_n,iyyy,month,iday)

hr=(date_n - idate_n)*24.0

print '(A,i5,i3,i3,2x,f3.0)','iyyy,month,iday,hr= ',&

& iyyy,month,iday,hr

is=1;ie=len("output/")

write(ofle(is:ie),'(A)')"output/"

is=ie+1;ie=is+lnblnk(varname)-1

write(ofle(is:ie),'(A)')varname(1:lnblnk(varname))

is=ie+1;ie=is+10

write(ofle(is:ie),'(i4.4,2i2.2 ,A, i2.2)')iyyy,month,iday,'_',int(hr)

! print '(A)',ofle(1:lnblnk(ofle))

is=ie+1;ie=is+len(".asc")-1

write(ofle(is:ie),'(A)')".asc"

print '(A)',ofle(1:lnblnk(ofle))

print '(A,A)','Output: ',ofle(1:lnblnk(ofle))

print *

open(20,file=ofle)

write(20,'(A,A)') '# input= ',infle(1:lnblnk(infle))

write(20,'(A,i4.4,2i2.2,i3.2,A)')'# yyyymmdd,hr= ',&

& iyyy,month,iday,int(hr),' GMT'

do j=1,jm

do i=1,im

varout=sdata(i,j,klevel,nstep)*scale + offset

write(20,'(2f12.6,f12.4,1x,2i5)') &

lon(i),lat(j),varout,i,j

enddo !i

enddo !j

close(20)

! *1： scaleとは、ある変数の平均的なオーダー

! *2： offsetとは、ある変数の平均的な値

! data(i,j) = sdata(i,j) * scale + offset

! のようにして、本当のデータを得る

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

write(*,*)

end program read_ncep_4d

----------------------

End of src/read_ncep_4d_int.f90

----------------------

======================

src/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=.

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

#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 ${OBJDIR}

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

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

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

LDFLAGS=-L/usr/local/netcdf-4.1.3/lib -lnetcdf -lnetcdff -I/usr/local/netcdf-4.1.3/include -module ${OBJDIR}

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

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

TARGET1=../read_ncep_4d

TARGET2=

TARGETS=$(TARGET1) $(TARGET2)

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

MOD=

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

SRC7=

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

SRC9=read_ncep_4d.f90 \

caldat.f90 julday.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

----------------------

End of src/makefile

----------------------