MGDSSTデータをMSMデータの格子に内挿する

実行用シェルスクリプト

#!/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

&para

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