OISSTデータから気候値を計算する

OISSTデータから指定された日の長期間平均を計算して結果をテキストファイルに書き出します。

出力範囲（緯度・経度）も指定できます。

下記の例では、1985年から2004年までの期間の７月31日における平均値を計算します。

実行の際のパラメータの設定は実行用シェルスクリプト(oisstclim.sh)の中で指定します。

コンパイル用のmakefileはifortでないと上手く動かないかもしれません。makefileを自分でかけばgfortranでも動くと思います。

実行環境

sh-3.2$ cat /etc/redhat-release

CentOS release 5.9 (Final)

sh-3.2$ ifort --version

ifort (IFORT) 13.1.1 20130313

sh-3.2$ GMT --version

GMT Version 4.4.0

Copyright 1991-2009 Paul Wessel and Walter H. F. Smith

ソースファイル

下記のファイルをsrcというサブディレクトリに置く。

sh-3.2$ ls src/*.f90 src/makefile

src/makefile src/mod_output.f90 src/scaling.f90

src/mod_ave_time.f90 src/mod_size_array.f90 src/setfillvalue.f90

src/mod_namelist.f90 src/mod_time_control.f90

src/mod_netcdfio.f90 src/oisstclim.f90

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

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 = #-convert big_endian -module ${OBJDIR}

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

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

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

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

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

TARGET1=../oisstclim

TARGET2=

TARGETS=$(TARGET1) $(TARGET2)

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

MOD=mod_namelist.f90 mod_size_array.f90 mod_netcdfio.f90 ¥

mod_time_control.f90 ¥

mod_output.f90

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

SRC7=

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

SRC9=oisstclim.f90 ¥

setfillvalue.f90 scaling.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

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

mod_ave_time.f90

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

module mod_ave_time

! Description:

!

! Author: am

!

! Host: aofd30

! Directory: /work2/am/WRF/TimeAve/src

!

! Revision history:

! This file is created by /usr/local/mybin/nff.sh at 16:07 on 02-10-2013.

! use

! implicit none

contains

subroutine ave_time_var2d(ncid,varid,it0,it1,ni,nj,var2d_tav)

integer,intent(in)::ncid, varid

integer,intent(in)::it0,it1

integer,intent(in)::ni,nj

real, intent(inout)::var2d_tav(ni,nj)

real, allocatable:: varin(:,:,:)

integer error

integer start(3)

integer dcount(3)

include 'netcdf.inc'

nt=it1-it0+1

allocate(varin(ni,nj,nt))

start(1) = 1

start(2) = 1

start(3) = it0

dcount(1) = ni

dcount(2) = nj

dcount(3) = nt

! dcount(3) = 1

error = nf_get_vara_real(ncid,varid,start,dcount,varin)

if (error /= 0) then

print *,'Error at nf_get_vara_real in ave_time_var2d (module mod_ave_time)'

print *,'ncid = ' ,ncid

print *,'varid = ' ,varid

print *,'start = ' ,start

print *,'dcount = ' ,dcount

stop

endif

do n=1,nt

var2d_tav(:,:)=var2d_tav(:,:)+varin(:,:,n)

enddo !n

! Here I used a stupid method because I want to avoid a mistake.

! do i=it0,it1

! start(3)=i

! error = nf_get_vara_real(ncid,varid,start,dcount,varin2d)

! if (error /= 0) then

! print *,'Error at nf_get_vara_real in ave_time_var2d (module mod_ave_time)'

! print *,'i = ' ,i

! print *,'ncid = ' ,ncid

! print *,'varid = ' ,varid

! print *,'start = ' ,start

! print *,'dcount = ' ,dcount

! stop

! endif

! var2d_tav(:,:)=var2d_tav(:,:)+varin2d(:,:)

! enddo !i

var2d_tav(:,:)=var2d_tav(:,:)/float(it1-it0+1)

end subroutine ave_time_var2d

end module mod_ave_time

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

mod_namelist.f90

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

module mod_namelist

! Description:

!

! Author: am

!

! Host: aofd30

! Directory: /work2/am/WRF/TimeAve/src

!

! Revision history:

! This file is created by /usr/local/mybin/nff.sh at 18:33 on 02-08-2013.

! use

! implicit none

! write(*,'(a)')'Module: mod_namelist'

! write(*,*)''

character(len=500),save::infle

character(len=500),save::ofle

character(len=4), save ::mmdd

integer,save::yyyy1,yyyy2

real,save:: lon1, lon2, lat1, lat2

contains

subroutine check_namelist

print '(a )','subroutine check_namelist'

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

print '(a,a )','ofle = ',trim(ofle)

print '(a,a )','mmdd= ',trim(mmdd)

print '(a,i5)','yyyy1= ',yyyy1

print '(a,i5)','yyyy2= ',yyyy2

print '(a,f9.3)','lon1= ',lon1

print '(a,f9.3)','lon2= ',lon2

print '(a,f9.3)','lat1= ',lat1

print '(a,f9.3)','lat2= ',lat2

print '(a )','end subroutine check_namelist'

print *

end subroutine check_namelist

! write(*,'(a)')'Done module mod_namelist.'

! write(*,*)

end module mod_namelist

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

mod_netcdfio.f90

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

module mod_netcdfio

! Description:

!

! Author: am

!

! Host: aofd30

! Directory: /work2/am/WRF/TimeAve

!

! Revision history:

! This file is created by /usr/local/mybin/nff.sh at 20:07 on 02-08-2013.

integer maxdim,nvars,ngatts,nunlimit

contains

subroutine read_nc_header(infle, ncid)

include 'netcdf.inc'

character(len=*),intent(in)::infle

integer,intent(out)::ncid ! netCDF ID number

integer error

character(len=80) dimname

! print '(A)','subroutine read_nc_header'

error = nf_open(infle,NF_NOWRITE,ncid)

if (error /= 0) then

print *,'Error at nf_open in read_nc_header (module mod_netcdfio)'

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

stop

endif

error = nf_inq(ncid,maxdim,nvars,ngatts,nunlimit)

if (error /= 0) then

print *,'Error at nf_inq in read_nc_header (module mod_netcdfio)'

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

print *,'ncid = ',ncid

stop

endif

! print *,'maxdim = ',maxdim

! print *,'nvars = ',nvars

! print *,'ngatts = ',ngatts

! print *,'nunlimit = ',nunlimit

! print '(A)','end subroutine read_nc_header'

! print *

end subroutine read_nc_header

subroutine read_nc_varid(ncid, varname, varid)

integer,intent(in)::ncid

character(len=*),intent(in)::varname

integer,intent(out)::varid

include 'netcdf.inc'

! print '(A)','read_nc_varid'

error = nf_inq_varid(ncid,varname,varid)

if (error /= 0) then

print *,'Error at nf_inq_varid in read_nc_varid (module mod_netcdfio)'

print *,'ncid = ',ncid

print *,'varname= ',varname

stop

endif

! print *,'varname= ',varname

! print *,'varid = ',varid

! print '(A)','end subroutine read_nc_varid'

! print *

end subroutine read_nc_varid

subroutine read_nc_size_array(ncid, varid, n)

integer,intent(in)::ncid, varid

integer,intent(out)::n

include 'netcdf.inc'

integer error

character(len=80) dimname

integer dimids(10)

integer i,j

integer maxdim,nvars,ngatts,nunlimit

integer ndims

! print '(A)','subroutine read_nc_size_array'

n=0

error = nf_inq_varndims(ncid,varid,ndims)

if (error /= 0) then

print *,'Error at nf_inq_varndims in read_nc_size_array (module mod_netcdfio)'

print *,'ncid = ',ncid

print *,'varid = ',varid

stop

endif

error = nf_inq_vardimid(ncid,varid,dimids)

if (error /= 0) then

print *,'Error at nf_inq_vardimid in read_nc_size_array (module mod_netcdfio)'

print *,'ncid = ',ncid

print *,'varid = ',varid

stop

endif

do i=1,ndims

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

if (error /= 0) then

print *,'Error at nf_inq_dimname in read_nc_size_array (module mod_netcdfio)'

print *,'ncid = ',ncid

print *,'i = ',i

print *,'dimids = ',dimids(i)

stop

endif

error = nf_inq_dimlen(ncid,dimids(i),n)

if (error /= 0) then

print *,'Error at nf_inq_dimlen in read_nc_size_array (module mod_netcdfio)'

print *,'ncid = ',ncid

print *,'dimname= ',dimname

print *,'i = ',i

print *,'dimids = ',dimids(i)

stop

endif

enddo

! print '(A)','end subroutine read_nc_size_array'

! print *

end subroutine read_nc_size_array

subroutine read_nc_float1d(ncid, varid, n, var1real)

integer,intent(in)::ncid, varid

integer,intent(in)::n

real,intent(inout)::var1real(n)

integer error

include 'netcdf.inc'

! print '(A)','subroutine read_nc_2dvar'

error = nf_get_var_real(ncid,varid,var1real)

if (error /= 0) then

print *,'Error at nf_get_var_real in read_nc_float1d (module mod_netcdfio)'

print *,'ncid = ' ,ncid

print *,'varid = ' ,varid

print *,'start = ' ,start

print *,'dcount = ' ,dcount

stop

endif

! print *,'ncid = ' ,ncid

! print *,'varid = ' ,varid

! print *,'start = ' ,start

! print *,'dcount = ' ,dcount

! print '(A)','end subroutine read_nc_float1d'

! print *

end subroutine read_nc_float1d

subroutine read_nc_short2d(ncid, varid, ni, nj, it, varin2d)

integer,intent(in)::ncid, varid

integer,intent(in)::ni,nj

integer,intent(in)::it

real,intent(inout)::varin2d(ni,nj)

integer(2)::var2dint2(ni,nj)

integer error

integer start(3)

integer dcount(3)

include 'netcdf.inc'

! print '(A)','subroutine read_nc_short2d'

start(1) = 1

start(2) = 1

start(3) = it

dcount(1) = ni

dcount(2) = nj

dcount(3) = 1

error = nf_get_vara_int2(ncid,varid,start,dcount,var2dint2)

if (error /= 0) then

print *,'Error at nf_get_vara_real in read_nc_short2d (module mod_netcdfio)'

print *,'ncid = ' ,ncid

print *,'varid = ' ,varid

print *,'start = ' ,start

print *,'dcount = ' ,dcount

stop

endif

! print *,'ncid = ' ,ncid

! print *,'varid = ' ,varid

! print *,'start = ' ,start

! print *,'dcount = ' ,dcount

varin2d(:,:)=float(var2dint2(:,:))

! print '(A)','end subroutine read_nc_short2d'

! print *

end subroutine read_nc_short2d

subroutine read_nc_scales(ncid,varid,scale,offset)

integer,intent(in)::ncid,varid

real,intent(inout)::scale,offset

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

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

end subroutine read_nc_scales

! write(*,'(a)')'Done module mod_netcdfio.'

! write(*,*)

end module mod_netcdfio

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

mod_output.f90

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

module mod_output

! Description:

!

! Author: am

!

! Host: aofd30

! Directory: /work2/am/WRF/TimeAve

!

! Revision history:

! This file is created by /usr/local/mybin/nff.sh at 13:26 on 02-10-2013.

! use

! implicit none

! write(*,'(a)')'Module: mod_output'

! write(*,*)''

contains

subroutine write2d(infle,ofle,mmdd,varname, &

ni,nj,lon,lat,lon1,lon2,lat1,lat2,var2d)

! Description:

! Author: am

! Host: aofd30

! Directory: /work2/am/WRF/TimeAve/src

! Revision history:

! This file is created by /usr/local/mybin/nff.sh at 12:25 on 02-10-2013.

character(len=*),intent(in)::infle,ofle,mmdd,varname

integer,intent(in)::ni,nj

real,intent(in)::lon1,lon2,lat1,lat2

real,intent(in)::lon(ni),lat(nj),var2d(ni,nj)

write(*,'(a)')'Subroutine: write2d'

! Set output file name

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

open(iu,file=ofle)

call print_header(iu)

call print_ioinfo(iu,infle,"input")

! write(iu,'(A,A)') '# varname= ',varname

write(iu,'(A,A)')'# mmdd= ', mmdd

! write(iu,'(A,A)')'# lon, lat, ',varname

do j=1,nj

if(lat(j) >= lat1)then

js=j

exit

endif

enddo !j

do j=nj,1,-1

if(lat(j) <= lat2)then

je=j

exit

endif

enddo !j

do i=1,ni

if(lon(i) >= lon1)then

is=i

exit

endif

enddo !j

do i=ni,1,-1

if(lon(i) <= lon2)then

ie=i

exit

endif

enddo !j

print *,is,ie,js,je

do j=js,je

do i=is,ie

write(iu,'(f12.6,f11.6,e14.6)')lon(i),lat(j),var2d(i,j)

enddo !i

enddo !j

close(iu)

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

! write(*,*)

end subroutine write2d

subroutine print_header(iunit)

integer,intent(in):: iunit

character(len=500) :: dirinfo, userinfo, hostinfo

INTEGER date_time(8)

CHARACTER REAL_CLOCK(3)*12

character(len=500) :: progname,shellinfo

! !REMARKS:

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

!

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

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

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

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

! 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: ',hostinfo(1:lnblnk(hostinfo))

write(iunit,'(A,A)')'# cwd: ',dirinfo(1:lnblnk(dirinfo))

write(iunit,'(A,A)')'# user: ',userinfo(1:lnblnk(userinfo))

! Get program name

call getarg( 0, progname)

write(iunit,'(A,A)')'# program name: ',progname(1:lnblnk(progname))

return

end subroutine print_header

subroutine print_ioinfo(iunit,fle,iokind)

integer,intent(in):: iunit

character,intent(in):: fle*(*), iokind*(*)

write(iunit,'(A,A,A,A)')'# ',iokind,': ',fle(1:lnblnk(fle))

end subroutine print_ioinfo

! write(*,'(a)')'Done module mod_output.'

! write(*,*)

end module mod_output

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

mod_size_array.f90

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

module mod_size_array

! Description:

!

! Author: am

!

! Host: aofd30

! Directory: /work2/am/WRF/TimeAve/src

!

! Revision history:

! This file is created by /usr/local/mybin/nff.sh at 22:00 on 02-08-2013.

! use

! implicit none

integer,save:: ni,nj,nt

end module mod_size_array

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

mod_time_control.f90

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

module mod_time_control

! Description:

!

! Author: am

!

! Host: aofd30

! Directory: /work2/am/WRF/TimeAve

!

! Revision history:

! This file is created by /usr/local/mybin/nff.sh at 13:36 on 02-10-2013.

! use

! implicit none

! write(*,'(a)')'Module: mod_time_control'

! write(*,*)''

contains

subroutine find_time_step(yyyy1,mmdd,it)

integer,intent(in)::yyyy1

character(len=*),intent(in)::mmdd

integer,intent(inout)::it

integer yyyy,mm,dd

integer yyyy0,mm0,dd0,yearday

read(mmdd(1:2),*)mm

read(mmdd(3:4),*)dd

yyyy0=yyyy1

mm0=1

dd0=1

yyyy=yyyy1

!print *,yyyy,mm,dd

!print *,yyyy0,mm0,dd0

yearday=julday(yyyy,mm,dd)-julday(yyyy0,mm0,dd0)+1

it=yearday

end subroutine find_time_step

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.

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

!

INTEGER IGREG

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

INTEGER JY,JM,JA

!

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 module mod_time_control

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

oisstclim.f90

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

program oisstclim

! Description:

!

! Author: am

!

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

! Directory: /work1/am/SST/OISST/Seasonal/src.test

!

! Revision history:

! This file is created by /usr/local/mybin/nff.sh at 17:14 on 04-24-2013.

use mod_namelist

use mod_size_array

use mod_netcdfio

use mod_time_control

!use mod_ave_time

use mod_output

character(len=500):: input_file

integer varidlon, varidlat, varidtime, varidsst

! Variable ID used in netCDF file

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

character(len=19)::time !e.g.,"2012-07-11_00:00:00"

real,allocatable::var2d_in(:,:) !For storing 2d input data

real,allocatable::ave2d(:,:) !Time average of 2d data data

real threshold, fillvalue

namelist/para/infle, ofle, mmdd, yyyy1, yyyy2, lon1, lon2, lat1, lat2, &

threshold,fillvalue

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

write(*,*)''

read(*,para)

call check_namelist

YEAR: do n=yyyy1,yyyy2

input_file=infle

is=index(infle,'????'); ie=is+3

write(input_file(is:ie),'(i4.4)')n

print *,trim(input_file)

! Get netCDF file ID

call read_nc_header(input_file, ncid)

! Read lon and lat

call read_nc_varid(ncid, 'lon' , varidlon)

call read_nc_varid(ncid, 'lat' , varidlat)

call read_nc_varid(ncid, 'time', varidtime)

call read_nc_varid(ncid, 'sst' , varidsst)

call read_nc_size_array(ncid, varidlon, ni)

! print *,'ni=',ni

call read_nc_size_array(ncid, varidlat, nj)

! print *,'nj=',nj

call read_nc_size_array(ncid, varidtime, nt)

! print *,'nt=',nt

if(n == yyyy1)then

allocate(lon(ni),lat(nj),var2d_in(ni,nj),ave2d(ni,nj))

ave2d(:,:)=0.0

endif

call read_nc_float1d(ncid, varidlon, ni, lon)

call read_nc_float1d(ncid, varidlat, nj, lat)

call read_nc_scales(ncid,varidsst,scale,offset)

! print *,'scale= ',scale

! print *,'offset=',offset

call find_time_step(yyyy1,mmdd,it)

! print *,'it=',it

call read_nc_short2d(ncid, varidsst, ni, nj, it, var2d_in)

ave2d(:,:)=ave2d(:,:)+var2d_in(:,:)

enddo YEAR

d=float(yyyy2-yyyy1+1)

print *,'d= ',d

ave2d(:,:)=ave2d(:,:)/d

call scaling(ni,nj,scale,offset,ave2d)

call setfillvalue(ni,nj,threshold,fillvalue,ave2d)

call write2d(input_file,ofle,mmdd,'sst',&

&ni,nj,lon,lat,lon1,lon2,lat1,lat2,ave2d)

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

write(*,*)

end program oisstclim

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

scaling.f90

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

subroutine scaling(ni,nj,scale,offset,var2d)

! Description:

!

! Author: am

!

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

! Directory: /work1/am/WRFpar/PreProc/OISST.Daily.Climat/src

!

! Revision history:

! This file is created by /usr/local/mybin/nff.sh at 18:56 on 06-04-2013.

! use

! implicit none

integer,intent(in)::ni,nj

real,intent(in) ::scale,offset

real,intent(inout)::var2d(ni,nj)

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

! write(*,*)''

var2d(:,:)=var2d(:,:)*scale+offset

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

! write(*,*)

end subroutine scaling

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

setfillvalue.f90

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

subroutine setfillvalue(ni,nj,threshold,fillvalue,var2d)

! Description:

!

! Author: am

!

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

! Directory: /work1/am/WRFpar/PreProc/OISST.Daily.Climat/src

!

! Revision history:

! This file is created by /usr/local/mybin/nff.sh at 18:44 on 06-04-2013.

! use

! implicit none

integer,intent(in)::ni,nj

real,intent(in) ::threshold,fillvalue

real,intent(inout)::var2d(ni,nj)

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

! write(*,*)''

do j=1,nj

do i=1,ni

if(var2d(i,j) > threshold)then

var2d(i,j)=fillvalue

endif

enddo

enddo

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

! write(*,*)

end subroutine setfillvalue

コンパイル

sh-3.2$ cd src

sh-3.2$ make clean; make

rm -rf ../obj/mod_namelist.o ../obj/mod_size_array.o ../obj/mod_netcdfio.o ../obj/mod_time_control.o ../obj/mod_output.o ../obj/oisstclim.o ../obj/setfillvalue.o ../obj/scaling.o ../obj/*.mod

if [ ! -d ../obj ]; then \

mkdir -p ../obj ; \

fi

ifort -c -O3 -module ../obj -c -o ../obj/mod_namelist.o mod_namelist.f90

ifort -c -O3 -module ../obj -c -o ../obj/mod_size_array.o mod_size_array.f90

ifort -c -O3 -module ../obj -c -o ../obj/mod_netcdfio.o mod_netcdfio.f90

ifort -c -O3 -module ../obj -c -o ../obj/mod_time_control.o mod_time_control.f90

ifort -c -O3 -module ../obj -c -o ../obj/mod_output.o mod_output.f90

ifort -c -O3 -module ../obj -c -o ../obj/oisstclim.o oisstclim.f90

ifort -c -O3 -module ../obj -c -o ../obj/setfillvalue.o setfillvalue.f90

ifort -c -O3 -module ../obj -c -o ../obj/scaling.o scaling.f90

ifort -o ../oisstclim ../obj/mod_namelist.o ../obj/mod_size_array.o ../obj/mod_netcdfio.o ../obj/mod_time_control.o ../obj/mod_output.o ../obj/oisstclim.o ../obj/setfillvalue.o ../obj/scaling.o -module ../obj -L/usr/local/lib -lnetcdf

実行用シェルスクリプト

sh-3.2$ cat oisstclim.run.sh

#!/bin/bash

# Description:

#

# Author: am

#

# Host: aofd30

# Directory: /work2/am/WRF/TimeAve

#

# Revision history:

# This file is created by /usr/local/mybin/nbscr.sh at 16:27 on 02-06-2013.

usage(){

cat << EOF

Usage: $0 [-d <indir>] [-i <input>] [-o <output>] [-t <mmdd>] \

[-f <yyyy1>] [-l <yyyy2>] \

[-e <lon1>] [-w <lon2>] [-s <lat1>] [-n <lat2>] \

[-l <threshold>] [-v <fillvalue>] \

[-h]

-d: Set input directory to <indir>

<indir> = input directory (e.g., input)

-i: Set input file name to <input>

<input> = input file name (e.g., sst.day.mean.1982.v2.nc)

-o: Set output file name to <output>

<output> = output file name (e.g., sst.day.mean.v2.19820101.txt)

-t: set date in 4 digits (e.g., Jan., 1st -> 0101 )

-f: First year to be used for computing climatology (e.g., 1982)

-l: Last year to be used for computing climatology (e.g., 2011)

-e: set longitude of east end of the output file to <lon1> [degree]

-w: set longitude of west end of the output file to <lon2> [degree]

-s: set latitude of south end of the output file to <lat1> [degree]

-n: set latitude of north end of the output file to <lat2> [degree]

-l: threshold value

-v: fillvalue

-h: Print this help

EOF

}

# Default values

flaglog="false"

flagindir="false"

indir="input"

input="sst.day.mean.????.v2.nc"

mmdd="0731"

yyyy1="1985"

yyyy2="2004"

outdir="output"

lon1=104 #100 #0

lon2=146 #180 #360

lat1=15 #0 #-90

lat2=50 #50 #90

threshold=327.

fillvalue=-0.999000E+09

while getopts d:i:o:t:f:l:w:e:s:n:l:v:h OPT; do

case $OPT in

"d" ) flagindir="true" ;indir=$OPTARG ;;

"i" ) flaginput="true" ;input=$OPTARG ;;

"o" ) flagoutput="true" ;output=$OPTARG ;;

"t" ) flagtime="true" ;mmdd=$OPTARG ;;

"f" ) flagyyyy1="true" ;yyyy1=$OPTARG ;;

"l" ) flagyyyy2="true" ;yyyy2=$OPTARG ;;

"w" ) flagwest="true" ;lon1=$OPTARG ;;

"e" ) flageast="true" ;lon2=$OPTARG ;;

"s" ) flagsouth="true" ;lat1=$OPTARG ;;

"n" ) flagnorth="true" ;lat2=$OPTARG ;;

"l" ) flagthold="true" ;threshold=$OPTARG ;;

"v" ) flagfillv="true" ;fillvalue=$OPTARG ;;

"h" ) usage ;exit 0 ;;

* ) usage

esac

done

shift $(expr $OPTIND - 1)

output="${outdir}/sst.clim.${mmdd}_${yyyy1}-${yyyy2}.txt"

echo

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

echo

exe=oisstclim

if [ ! -f $exe ]; then

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

exit 1

fi

if [ ! -d $outdir ]; then

echo Creating directory, $outdir

mkdir -p $outdir

fi

echo -------------------------------------------------------------------------

echo

echo Set-up input file name

echo

if [ ! -d $indir ]; then

echo No such directory, $indir

exit 1

fi

infle=${indir}/${input}

if [ ! -f $infle ]; then

echo No such file, $infle

exit 1

fi

echo -------------------------------------------------------------------------

echo

echo Set-up namelist file

echo

namelist=namelist.$(basename $0 .sh).txt

cat << EOF > $namelist

&para

infle="${infle}"

ofle="${output}"

mmdd="${mmdd}"

yyyy1=${yyyy1}

yyyy2=${yyyy2}

lon1=${lon1}

lon2=${lon2}

lat1=${lat1}

lat2=${lat2}

threshold=${threshold}

fillvalue=${fillvalue}

&end

EOF

echo ${namelist}:

cat $namelist

echo

echo Done.

echo -------------------------------------------------------------------------

echo -------------------------------------------------------------------------

echo

echo Run ${exe}

echo

${exe} < $namelist

echo

echo "Done ${exe}"

echo

echo -------------------------------------------------------------------------

echo

echo "Done $0"

echo

実行例

sh-3.2$ oisstclim.run.sh

Shell script, oisstclim.run.sh starts.

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

Set-up input file name

./oisstclim.run.sh: line 111: [: too many arguments

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

Set-up namelist file

namelist.oisstclim.run.txt:

&para

infle="input/sst.day.mean.????.v2.nc"

ofle="output/sst.clim.0731_1985-2004.txt"

mmdd="0731"

yyyy1=1985

yyyy2=2004

lon1=104

lon2=146

lat1=15

lat2=50

threshold=327.

fillvalue=-0.999000E+09

&end

Done.

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

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

Run oisstclim

Program oisstclim starts.

subroutine check_namelist

infle = input/sst.day.mean.????.v2.nc

ofle = output/sst.clim.0731_1985-2004.txt

mmdd= 0731

yyyy1= 1985

yyyy2= 2004

lon1= 104.000

lon2= 146.000

lat1= 15.000

lat2= 50.000

end subroutine check_namelist

input/sst.day.mean.1985.v2.nc

input/sst.day.mean.1986.v2.nc

input/sst.day.mean.1987.v2.nc

input/sst.day.mean.1988.v2.nc

input/sst.day.mean.1989.v2.nc

input/sst.day.mean.1990.v2.nc

input/sst.day.mean.1991.v2.nc

input/sst.day.mean.1992.v2.nc

input/sst.day.mean.1993.v2.nc

input/sst.day.mean.1994.v2.nc

input/sst.day.mean.1995.v2.nc

input/sst.day.mean.1996.v2.nc

input/sst.day.mean.1997.v2.nc

input/sst.day.mean.1998.v2.nc

input/sst.day.mean.1999.v2.nc

input/sst.day.mean.2000.v2.nc

input/sst.day.mean.2001.v2.nc

input/sst.day.mean.2002.v2.nc

input/sst.day.mean.2003.v2.nc

input/sst.day.mean.2004.v2.nc

d= 20.00000

Subroutine: write2d

Output: output/sst.clim.0731_1985-2004.txt

417 584 421 560

Done subroutine write2d.

Done program oisstclim.

Done oisstclim

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

Done ./oisstclim.run.sh

描画

sh-3.2$ cd plot

sh-3.2$ cat pl_sst.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

flagn=${flagn:-"NOT_AVAILABLE"}

shift `expr $OPTIND - 1`

. ./gmtpar.sh

export LANG=C

echo "Bash script $0 starts."

scale=1.

range1=116/140/23/42

size=4 #2.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="20/30/.5"

makecpt -Cpanoply -T$cpt_range > $cpt

plcontour=n

CI1=0.5

CI2=1

scalex=2

scaley=-0.5

scalesize=2.5

scale_anot=a5f1/:"${deg}C":

indir="../output/"

if [ $# -ne 1 ]; then

usage

exit 1

fi

in=$indir$1

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 .asc).ps

echo Input : $in

echo Output : $out

grd=$1.grd #$(basename $in .asc).grd

tstamp1=$(ls -l $in |awk '{print $6,$7,$8}')

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 -JQ135/$size \

-K -X1.5 -Y7 -P > $out

grdimage $grd -R$range1 -JQ \

-C$cpt \

-O -K \

>> $out

cl=0

cu=40

grdcontour $grd -R$range1 -J${map} -C${CI1} -W4/${white} \

-L${cl}/${cu} -O -K >> $out

grdcontour $grd -R$range1 -J${map} -A${CI2}f10 -W8/${white} \

-L${cl}/${cu} -O -K >> $out

psmask -C -O -K >> $out

if [ $flagn = "true" ]; then

pscoast -R$range1 -JQ -B${anot}WSne -Dl -W2 -O -K >> $out

else

pscoast -R$range2 -JQ -B${anot}WSne -Dl -W2 -G200 -O -K >> $out

fi

psscale -D$scalex/$scaley/$scalesize/0.1h -B${scale_anot} -E \

-C$cpt -O -K >> $out

range2=120/132/23/36

grdimage $grd -R$range2 -JQ125/$size \

-C$cpt \

-O -K -Y-5.5 \

>> $out

cl=0

cu=40

grdcontour $grd -R$range2 -J${map} -C${CI1} -W4/${white} \

-L${cl}/${cu} -O -K >> $out

grdcontour $grd -R$range2 -J${map} -A${CI2}f10 -W8/${white} \

-L${cl}/${cu} -O -K >> $out

#psmask -C -O -K >> $out

if [ $flagn = "true" ]; then

pscoast -R$range2 -JQ -B${anot}WSne -Dl -W2 -O -K >> $out

else

pscoast -R$range2 -JQ -B${anot}WSne -Dl -W2 -G200 -O -K >> $out

fi

psscale -D$scalex/$scaley/$scalesize/0.1h -B${scale_anot} -E \

-C$cpt -O -K >> $out

xoffset=-1

yoffset=8.5

comment=

. ./note.sh

rm -rf ${grd} ${cpt}

echo "Done $0.sh"

sh-3.2$ cat gmtpar.sh

gmtset MEASURE_UNIT INCH

gmtset LABEL_FONT_SIZE 18

gmtset HEADER_FONT_SIZE 20

gmtset ANOT_FONT_SIZE 18

gmtset TICK_PEN 4

# ?n?}?̏c?????ɎȁX???g??Ȃ?

gmtset BASEMAP_TYPE PLAIN

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

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"

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

gmtset D_FORMAT %lg #?f?t?H???g(?K??l)

#gmtset D_FORMAT %12.6f

#gmtset D_FORMAT %12.5g

sh-3.2$ cat 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}

sh-3.2$ pl_sst.sh sst.clim.0731_1985-2004.txt

Bash script ./pl_sst.sh starts.

Input : ../output/sst.clim.0731_1985-2004.txt

Output : ../ps/sst.clim.0731_1985-2004.txt.ps

pstext: Record 6 is incomplete (skipped)

Done ./pl_sst.sh.sh