実行用シェルスクリプト
#!/bin/sh
exe=interpolate_mgdsst
parafile=${exe}.para.txt
outdir=output
if [ ! -d $outdir ]; then
mkdir -p $outdir
fi
yyyy=2012
mm=06
dd=14
indir=/work4/DATA/SST/MGDSST/data/${yyyy}
in=mgd_sst_glb_D${yyyy}${mm}${dd}.txt
infle=$indir/$in
msmim=481
msmjm=505
cat << EOF > $parafile
¶
infle="$infle",
msmim=${msmim},
msmjm=${msmjm},
&end
EOF
$exe < $parafile
makefile
#
#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=interpolate_mgdsst
TARGET2=
TARGETS=$(TARGET1) $(TARGET2)
#j MOD (Fortran90以降のモジュールファイル)
MOD=
#j SRC7 (Fortran77のソースファイル)
SRC7=
#j SRC9 (Fortran90のソースファイル)
SRC9=interpolate_mgdsst.f90 basename.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
Fortranプログラム
program interpolate_mgdsst
! implicit none
character infle*500,ofle*500, ofle2*500, ofle2_base*500, ofle3*500
character*100 dimname
integer :: iyyyy,imm,idd
integer :: i,j,jj
integer,parameter :: im=1440,jm=720
integer :: isst(im,jm),ist
real :: sst(im,jm)
real :: undef = -1.e+35
integer :: siz,wlength
real,allocatable :: msmsst(:,:), msmlon(:),msmlat(:)
real,parameter:: lonw=0.0, lats=-89.875 ! Origin of MGDSST coordinates
real,parameter:: dx=0.25, dy=0.25 ! Grid spacings of MGDSST
real xa(4),ya(4),za(4)
namelist /para/infle,msmim,msmjm
read(*,nml=para)
allocate(msmsst(msmim,msmjm), msmlon(msmim),msmlat(msmjm))
print '(A,A)','Input: ',infle(1:lnblnk(infle))
open(55,file=infle,status='old',form='formatted',action="read")
read(55,*) iyyyy,imm,idd
write(6,*) iyyyy,imm,idd
do j=1,jm
read(55,10) (isst(i,j),i=1,im)
enddo
close(55)
10 format(1440i3)
do j=1,jm
do i=1,im
jj = jm -j +1
ist = isst(i,jj)
if(ist .eq. 888 .or. ist.eq.999) then
sst(i,j) = undef
else
sst(i,j) = real(ist)*0.1
endif
enddo
enddo
! Set grid systems in MSM
! allocate(msmsst(msmim,msmjm), msmlon(msmim), msmlat(msmjm) )
do j=1,msmjm
msmlat(j)= 47.6 - 0.05*float(j-1)
enddo !j
do i=1,msmim
msmlon(i)=120.0 + 0.0625*float(i-1)
enddo !i
! do j=1,msmjm
! print *,msmlat(j)
! enddo !j
! do i=1,msmim
! print *,msmlon(i)
! enddo !j
! stop
do j=1,msmjm
jsw=int((msmlat(j)-lats)/dy)+1
do i=1,msmim
isw=int((msmlon(i)-lonw)/dx)
! print *,isw,jsw, msmlon(i),msmlat(j)
xa(1)=lonw + dx*float(isw)
ya(1)=lats + dy*float(jsw)
za(1)=sst(isw,jsw)
xa(2)=xa(1)+dx
ya(2)=ya(1)
za(2)=sst(isw+1,jsw)
xa(3)=xa(2)
ya(3)=ya(2)+dy
za(3)=sst(isw+1,jsw+1)
xa(4)=xa(1)
ya(4)=ya(3)
za(4)=sst(isw,jsw+1)
if(za(1)==undef .or. za(2)==undef .or. za(3)==undef .or. &
za(4)==undef)then
msmsst(i,j)=undef
else
call bilnr(xa,ya,za, msmlon(i), msmlat(j), msmsst(i,j))
endif
! print *,msmlon(i), msmlat(j)
! print *,xa(1),ya(1),za(1)
! print *,xa(2),ya(2),za(2)
! print *,xa(3),ya(3),za(3)
! print *,xa(4),ya(4),za(4)
! print *,msmsst(i,j)
! print *
enddo !i
enddo !j
! inquire(iolength=wlength) sst(1,1)
! siz = jm * im * wlength
inquire(iolength=wlength) msmsst(1,1)
siz = msmjm * msmim * wlength
ofle2=""
is=1;ie=len('output/mgdsst.msm.')
write(ofle2(is:ie),'(A)')'output/mgdsst.msm.'
is=ie+1; ie=is+3
write(ofle2(is:ie),'(i4.4)')iyyyy
is=ie+1; ie=is+1
write(ofle2(is:ie),'(i2.2)')imm
is=ie+1; ie=is+1
write(ofle2(is:ie),'(i2.2)')idd
is=ie+1; ie=is+len('.bin')
write(ofle2(is:ie),'(a4)')'.bin'
print "(A,A)",ofle2(1:lnblnk(ofle2))
open(77,file=ofle2,access='direct',form='unformatted', &
& & status='replace',recl=siz)
write(77,rec=1) ((msmsst(i,j),i=1,msmim),j=1,msmjm)
close(77)
! print *,msmsst
if(imm == 1)month='Jan'
if(imm == 2)month='Feb'
if(imm == 3)month='Mar'
if(imm == 4)month='Apr'
if(imm == 5)month='May'
if(imm == 6)month='Jun'
if(imm == 7)month='Jul'
if(imm == 8)month='Aug'
if(imm == 9)month='Sep'
if(imm == 10)month='Oct'
if(imm == 11)month='Nov'
if(imm == 12)month='Dec'
is=1;ie=len('output/mgdsst.msm.')
write(ofle3(is:ie),'(A)')'output/mgdsst.msm.'
is=ie+1;ie=is+3
write(ofle3(is:ie),'(i4.4)')iyyyy
is=ie+1;ie=is+1
write(ofle3(is:ie),'(i2.2)')imm
is=ie+1;ie=is+len('.ctl')
write(ofle3(is:ie),'(a4)')'.ctl'
print '(A,A)','output: ',trim(ofle3)
print *
open(111,file=ofle3)
call basename(ofle2, ofle2_base, idxbi)
write(111,*) "dset ^mgdsst.msm.%y4%m2%d2.bin"
write(111,*) "title mgdsst"
write(111,*) "undef ",undef
write(111,*) "options template little_endian yrev"
write(111,*) "xdef ",msmim," linear 120.0 0.0625"
write(111,*) "ydef ",msmjm," linear 22.4 0.0500"
write(111,*) "zdef 1 linear 0.0 1"
write(111,'(A,A3,i4.4,A,A)') &
& "tdef 31 linear 00z01",month,iyyyy,' ','1dy'
write(111,*) "vars 1"
write(111,*) "sst 1 0 sst"
write(111,*) "endvars"
close(111)
write(*,'(a)')'Done program interpolate_mgdsst.'
write(*,*)
end program interpolate_mgdsst
!***********************************************************************
! bilnr.f
! TASK
! 2D-bilinear interpolation
! Reference
! http://www.library.cornell.edu/nr/bookfpdf/f3-6.pdf
!***********************************************************************
subroutine bilnr(xa,ya, za, x, y, z)
!----------------- PARAMETERS AND COMMON VARIABLES ---------------------
! implicit double precision(a-h,o-z)
!--------------------------- ARGUMENTS ---------------------------------
real, intent(in) :: xa(4), ya(4), za(4)
real, intent(in) :: x,y
real, intent(out) :: z
!------------------------- LOCAL VARIABLES -----------------------------
real t, u
!-----------------------------------------------------------------------
! Subscript, a indicates the input.
!
!
! For given xa(i),ya(i),za(i),x & y, z is estimated.
!
! (xa(4),ya(4),za(4)) (xa(3),ya(3),za(3))
! +------------------+
! | |
! | x |
! | (x,y,z) |
! | |
! | |
! | |
! | |
! | |
! +------------------+
! (xa(1),ya(1),za(1)) (xa(2),ya(2),za(2))
!
t=(x-xa(1))/(xa(2)-xa(1))
u=(y-ya(1))/(ya(4)-ya(1))
z = (1.-t)*(1.-u)*za(1) &
& + t*(1.-u)*za(2) &
& + t*u*za(3) &
& + (1.-t)*u*za(4)
return
end
!***********************************************************************
! .f
! NOTE
! Reference
!***********************************************************************
subroutine basename(infle, bifle, idxbi)
!------------------------- DECLARATION ---------------------------------
! implicit none
! implicit double precision(a-h,o-z)
!------------------------- INCLUDE FILE --------------------------------
! include 'comblk.h'
!-------------------------- PARAMETERS ---------------------------------
! parameter()
!------------------------- COMMON VARIABLES ----------------------------
! common //
!----------------------------ARGUMENTS ---------------------------------
! double precision
! real
! integer
! dimension
! Input
character infle*(*)
! Output
character bifle*(*)
integer idxbi
! Input/output
!------------------------- LOCAL VARIABLES -----------------------------
character testchar*1
!-----------------------------------------------------------------------
!---+$|--1----+----2----+----3----+----4----+----5----+----6----+----7--
last=index(infle,' ')-1
do i=last, 1, -1
testchar=infle(i:i)
id=index(testchar,'/')
if(id.ne.0)then
idx=i+1
goto 1000
end if
end do !i
! If there is no slash in infle
bifle=infle(1:last)
idxbi=last
! write(*,'(A)')infle
! write(*,'(A)')bifle
! write(*,*)idxbi
return
! If there is a slash in infle
1000 bifle=infle(idx:last)
idxbi=last-idx+1
! write(*,'(A)')infle
! write(*,'(A)')bifle
! write(*,*)idxbi
return
!---+$|--1----+----2----+----3----+----4----+----5----+----6----+----7--
end
作図用GrADSコントロールファイル
dset ^mgdsst.msm.%y4%m2%d2.bin
title mgdsst
undef -1.0000000E+35
options template little_endian yrev
xdef 481 linear 120.0 0.0625
ydef 505 linear 22.4 0.0500
zdef 1 linear 0.0 1
tdef 31 linear 00z01Jun2012 1dy
vars 1
sst 1 0 sst
endvars