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
¶
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:
¶
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