実行用シェルスクリプト: mgdsst.alongtrack.run.sh
#!/bin/bash
# Description:
#
# Host: aofd165.fish.nagasaki-u.ac.jp
# Directory: /work2/kunoki/to_manda_sensei/Tools/MGDSST.AlongTrack
#
# Revision history:
# This file is created by /usr/local/mybin/nbscr.sh at 17:28 on 09-11-2014.
exe=mgdsst.alongtrack
namelist=${exe}.namelist.txt
infle=MandA2012Leg3.track.txt
ofle=MGDSST.MandA2012Leg3.track.txt
mgd_dir=/work4/DATA/SST/MGDSST/data/2012/
mgd_prefix=mgd_sst_glb_D
mgd_postfix=txt
nlon=1440
nlat=720
slon=0.125
slat=-89.875
dlon=0.25
dlat=0.25
undef=-999.
cat <<EOF >$namelist
¶
infle="$infle",
ofle="$ofle",
mgd_dir="$mgd_dir", !/work4/DATA/SST/MGDSST/data/2012/
mgd_prefix="$mgd_prefix", !mgd_sst_glb_D
mgd_postfix="$mgd_postfix", !txt
undef=$undef
&end
&grid
nlon=$nlon
nlat=$nlat
slon=$slon
slat=$slat
dlon=$dlon
dlat=$dlat
&end
EOF
${exe} < $namelist
echo "Done $(basename $0)"
echo
船の航路データ(サンプル): MandA2012Leg3.track.txt
# date(yyyy,mm,dd) time_JST(hh,mm) lat(dd,mm.mmm) lon(ddd,mm.mmm) wh(m) wp(s) wd(deg) ws(knot) ta(degC) ts(degC) p(hPa) gyr ss
10 # total number of data
2012 06 16 23 49 00 26 31.145 126 51.953 1.3 10.0 153.000 14.3 27.60 27.5 1011.3 233.8 13.1
2012 06 16 23 50 00 26 31.020 126 51.751 1.3 10.0 152.500 14.3 27.70 27.5 1011.1 233.9 13.1
2012 06 16 23 51 00 26 30.895 126 51.549 1.3 10.0 150.600 14.8 27.70 27.5 1011.2 234.0 13.2
2012 06 16 23 52 00 26 30.771 126 51.347 1.3 10.0 150.800 14.4 27.70 27.5 1011.3 234.2 13.1
2012 06 16 23 53 00 26 30.648 126 51.144 1.2 10.0 149.300 14.6 27.70 27.4 1011.0 233.4 13.0
2012 06 16 23 54 00 26 30.523 126 50.943 1.2 10.0 150.700 14.4 27.70 27.5 1011.1 233.4 13.1
2012 06 16 23 55 00 26 30.398 126 50.743 1.2 10.0 150.000 14.7 27.70 27.5 1011.0 233.9 13.0
2012 06 16 23 56 00 26 30.275 126 50.544 1.2 10.0 151.200 14.6 27.70 27.5 1011.1 233.9 13.0
2012 06 16 23 57 00 26 30.152 126 50.343 1.1 10.0 150.300 14.3 27.70 27.5 1010.9 233.6 13.0
2012 06 16 23 58 00 26 30.029 126 50.144 1.2 10.0 149.700 14.4 27.70 27.5 1011.0 234.4 13.0
入力データの時刻はJSTで与える
UTCに変換された時刻が出力される
フォートランプログラム用makefile
#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}
#j ターゲット名(最終的に作りたい実行ファイル名)
TARGET1=mgdsst.alongtrack
TARGET2=
TARGETS=$(TARGET1) $(TARGET2)
#j MOD (Fortran90以降のモジュールファイル)
MOD=
#j SRC7 (Fortran77のソースファイル)
SRC7=
#j SRC9 (Fortran90のソースファイル)
SRC9=mgdsst.alongtrack.f90 \
mgdsst_lonlat.f90 read_mgdsst.f90 \
interpolate_space.f90 bilinear.f90 \
interpolate_time.f90 \
lst2utc.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
mgdsst.alongtrack.f90
program mgdsst_alongtrack
! Description:
!
! Host: aofd165.fish.nagasaki-u.ac.jp
! Directory: /work2/kunoki/to_manda_sensei/Tools/MGDSST.AlongTrack
!
! Revision history:
! This file is created by /usr/local/mybin/nff.sh at 17:11 on 09-11-2014.
! use
! implicit none
character infle*500,ofle*500
integer,allocatable,dimension(:)::yyyy,mm,dd,hh,mi
integer::ss
integer,allocatable,dimension(:)::uyyyy,umm,udd,uhh,umi
integer,parameter::tdh=+9
integer dd1,dd2
real,allocatable,dimension(:)::lond,lonm,latd,latm
real lon,lat
real,allocatable,dimension(:)::mlon,mlat
real,allocatable,dimension(:,:)::msst1,msst2
character(len=100)::mgd_dir,mgd_prefix,mgd_postfix
character(len=300)::mgd_file1,mgd_file2
!DEBUG
character(len=100)::out_check
integer :: siz,wlength
namelist /para/infle,ofle,mgd_dir,mgd_prefix,mgd_postfix,undef
namelist /grid/nlon,nlat,slon,slat,dlon,dlat
read(*,nml=para)
read(*,nml=grid)
print '(A,A)','Input : ',trim(infle)
print '(A,A)','Output: ',trim(ofle)
allocate(mlon(nlon),mlat(nlat))
call mgdsst_lonlat(nlon,nlat,slon,slat,dlon,dlat,mlon,mlat)
print *,mlon(1),mlon(nlon)
print *,mlat(1),mlat(nlat)
allocate(msst1(nlon,nlat),msst2(nlon,nlat))
open(10,file=infle,action="read")
read(10,*)
read(10,*)nt
allocate(yyyy(nt),mm(nt),dd(nt),hh(nt),mi(nt))
allocate(uyyyy(nt),umm(nt),udd(nt),uhh(nt),umi(nt))
allocate(lond(nt),lonm(nt),latd(nt),latm(nt))
do n=1,nt
read(10,*)yyyy(n),mm(n),dd(n),hh(n),mi(n),ss,latd(n),latm(n),&
& lond(n),lonm(n)
enddo !n
open(20,file=ofle)
do n=1,nt
mgd_file1=trim(mgd_dir)//trim(mgd_prefix)//"yyyymmdd"//"."//&
& trim(mgd_postfix)
mgd_file2=trim(mgd_dir)//trim(mgd_prefix)//"yyyymmdd"//"."//&
& trim(mgd_postfix)
call lst2utc(yyyy(n),mm(n),dd(n),hh(n), &
& uyyyy(n),umm(n),udd(n),uhh(n), tdh)
umi(n)=mi(n)
if(uhh(n) <= 12)then !Assuming SST is defined at 1200 UTC
dd1=udd(n)-1
else
dd1=udd(n)
endif
is=index(mgd_file1,"yyyymmdd")
write(mgd_file1(is:is+7),'(i4.4,i2.2,i2.2)')uyyyy(n),umm(n),dd1
dd2=dd1+1
is=index(mgd_file2,"yyyymmdd")
write(mgd_file2(is:is+7),'(i4.4,i2.2,i2.2)')uyyyy(n),umm(n),dd2
call read_mgdsst(mgd_file1,msst1,nlon,nlat,undef)
! ! DEBUG: Print SST
! write(out_check(1:8+8),'(i4.4,i2.2,i2.2,A)')uyyyy(n),umm(n),udd(n),'.sst.bin'
! inquire(iolength=wlength) msst1(1,1)
! siz = nlat * nlon * wlength
!
! open(20,file=out_check,access='direct',form='unformatted', &
! & status='replace',recl=siz)
! write(20,rec=1) ((msst1(i,j),i=1,nlon),j=1,nlat)
! close(20)
call read_mgdsst(mgd_file2,msst2,nlon,nlat,undef)
lon=lond(n)+lonm(n)/60.
lat=latd(n)+latm(n)/60.
call interpolate_space(msst1,nlon,nlat,mlon,mlat,lon,lat,sst1,undef)
call interpolate_space(msst2,nlon,nlat,mlon,mlat,lon,lat,sst2,under)
if(uhh(n) >= 12)then
time_s=(uhh(n)-12)*3600.0+umi(n)*60.0
else if(uhh(n) < 12)then
time_s=(uhh(n)+12)*3600.0+umi(n)*60.0
endif
call interpolate_time(sst1,sst2,time_s,sst,w1,w2,undef)
! DEBUG
! print *,'JST ',yyyy(n),mm(n),dd(n),hh(n)
! print *,'UTC ',uyyyy(n),umm(n),udd(n),uhh(n)
! print *,trim(mgd_file1)
! print *,trim(mgd_file2)
! print *,"interpolate_time: "
! print *,sst1,w1
! print *,sst2,w2
! print *,sst
! print *
write(20,'(i4.4,1x,i2.2,1x,i2.2,1x,i2.2,1x,i2.2, &
& f5.0,f8.4, f5.0,f8.4, 3f10.3)') &
& uyyyy(n),umm(n),udd(n),uhh(n),umi(n),&
& latd(n),latm(n),lond(n),lonm(n), &
& sst,sst1,sst2
! DEBUG
write(*,'(i4.4,1x,i2.2,1x,i2.2,1x,i2.2,1x,i2.2, &
& f5.0,f8.4, f5.0,f8.4, 3f10.3)') &
& uyyyy(n),umm(n),udd(n),uhh(n),umi(n),&
& latd(n),latm(n),lond(n),lonm(n), &
& sst,sst1,sst2
end do !n
close(20)
! write(*,'(a)')'Done program mgdsst.alongtrack.'
! write(*,*)
end program mgdsst_alongtrack
mgdsst_lonlat.f90
subroutine mgdsst_lonlat(nlon,nlat,slon,slat,dlon,dlat,mlon,mlat)
! Description:
!
! Author: kunoki
!
! Host: aofd165.fish.nagasaki-u.ac.jp
! Directory: /work2/kunoki/to_manda_sensei/Tools/MGDSST.AlongTrack
!
! Revision history:
! This file is created by /usr/local/mybin/nff.sh at 18:41 on 09-11-2014.
! use
! implicit none
integer,intent(in)::nlon,nlat
real,intent(in)::slon,slat,dlon,dlat
real,intent(inout)::mlon(nlon),mlat(nlat)
do i=1,nlon
mlon(i)=slon + dlon*float(i-1)
enddo
do i=1,nlat
mlat(i)=slat + dlat*float(i-1)
enddo
end subroutine mgdsst_lonlat
read_mgdsst.f90
subroutine read_mgdsst(mgd_file,msst,nlon,nlat,undef)
! Description:
!
! Author: kunoki
!
! Host: aofd165.fish.nagasaki-u.ac.jp
! Directory: /work2/kunoki/to_manda_sensei/Tools/MGDSST.AlongTrack
!
! Revision history:
! This file is created by /usr/local/mybin/nff.sh at 19:01 on 09-11-2014.
! use
! implicit none
character(len=*),intent(in)::mgd_file
real,intent(in)::undef
real,intent(inout)::msst(nlon,nlat)
integer::isst(nlon,nlat)
! print '(A)',trim(mgd_file)
open(11,file=mgd_file,action="read")
read(11,*)iyyy,imm,idd
! print *,iyyy,imm,idd
do j=nlat,1,-1
read(11,'(1440i3)') (isst(i,j),i=1,nlon)
do i=1,nlon
ist=isst(i,j)
if(ist .eq. 888 .or. ist.eq.999) then
msst(i,j) = undef
else
msst(i,j) = real(ist)*0.1 + 273.15
endif
enddo !i
enddo !j
close(11)
end subroutine read_mgdsst
interpolate_space.f90
subroutine interpolate_space(msst,nlon,nlat,mlon,mlat,lon,lat,sst,undef)
! Description:
!
! Author: kunoki
!
! Host: aofd165.fish.nagasaki-u.ac.jp
! Directory: /work2/kunoki/to_manda_sensei/Tools/MSM.profiles.at.Stations/src
!
! Revision history:
! This file is created by /usr/local/mybin/nff.sh at 11:34 on 06-03-2014.
implicit none
integer,intent(in)::nlon,nlat
real,intent(in)::msst(nlon,nlat)
real,intent(in)::mlon(nlon),mlat(nlat)
real,intent(in)::lon,lat
real,intent(inout)::sst
real,intent(in)::undef
real dlon,dlat
integer idx,jdx,ierr
integer n
real xa(4), ya(4), za(4)
real x,y
real z
integer::irr=0
dlon=mlon(2)-mlon(1) !Assuming uniform grid spacing
dlat=mlat(2)-mlat(1)
idx=(lon - mlon(1))/dlon + 1
jdx=(lat - mlat(1))/dlat + 1
xa(1)=mlon(idx); ya(1)=mlat(jdx+1)
xa(2)=mlon(idx+1); ya(2)=mlat(jdx+1)
xa(3)=mlon(idx+1); ya(3)=mlat(jdx)
xa(4)=mlon(idx); ya(4)=mlat(jdx)
za(1)=msst(idx, jdx+1)
za(2)=msst(idx+1, jdx+1)
za(3)=msst(idx+1, jdx )
za(4)=msst(idx, jdx )
ierr=0
do n=1,4
if(za(n) == undef)then
ierr=1
endif
enddo
if(ierr==1)then
sst=undef
return
endif
call bilinear(xa,ya, za, lon, lat, z)
sst=z
! print *,"interpolate_space: "
! print *,"1 ",xa(1),ya(1),za(1)
! print *,"2 ",xa(2),ya(2),za(2)
! print *,"3 ",xa(3),ya(3),za(3)
! print *,"4 ",xa(4),ya(4),za(4)
! print *,"A ",lon,lat,z
! print *
end subroutine interpolate_space
bilinear.f90
!***********************************************************************
! bilinear.f90
! TASK
! 2D-bilinear interpolation
! http://imagingsolution.blog107.fc2.com/blog-entry-142.html
! http://www.library.cornell.edu/nr/bookfpdf/f3-6.pdf
!***********************************************************************
subroutine bilinear(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
interpolate_time.f90
subroutine interpolate_time(x1,x2,t,x,w1,w2,undef)
! Description:
!
! Author: kunoki
!
! Host: aofd165.fish.nagasaki-u.ac.jp
! Directory: /work2/kunoki/to_manda_sensei/Tools/MSM.profiles.at.Stations/src
!
! Revision history:
! This file is created by /usr/local/mybin/nff.sh at 14:14 on 06-03-2014.
real,intent(in)::x1,x2,t
real,intent(in)::undef
real,intent(inout)::x,w1,w2
if(x1 == undef .or. x2 == undef)then
x=undef
return
endif
w1=( t )/86400.
w2=( 86400. - t )/86400.
x =x1*w1 + x2*w2
end subroutine interpolate_time
lst2utc.f90
subroutine lst2utc(lyr,lmo,ldy, lhr, uyr, umo, udy, uhr, tdh)
! Description:
!
! Author: am
!
! Host: aofd165.fish.nagasaki-u.ac.jp
! Directory: /work2/kunoki/to_manda_sensei/Tools/MGDSST.AlongTrack
!
! Revision history:
! This file is created by /usr/local/mybin/nff.sh at 09:19 on 09-12-2014.
! Reference
! 1行で書けるうるう年判別法
! http://d.hatena.ne.jp/Kappuccino/20081025/1224906495
!
implicit none
integer,intent(in):: lyr,lmo,ldy, lhr
integer,intent(inout)::uyr,umo,udy, uhr
integer,intent(in):: tdh
integer y
integer,dimension(12)::last_day
data last_day/31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31/
uyr=lyr
umo=lmo
udy=ldy
uhr=lhr
uhr=lhr-tdh
if(uhr<0)then
uhr=uhr+24
udy=udy-1
endif
if(udy<=0)then
umo=umo-1
udy=last_day(umo)
if(umo == 2)then
y=uyr
udy = 28 + (1 / (y - y / 4 * 4 + 1)) * &
& (1 - 1 / (y - y / 100 * 100 + 1)) &
& + (1 / (y - y / 400 * 400 + 1));
endif
endif
if(umo<=0)then
umo=1
uyr=uyr-1
endif
end subroutine lst2utc
出力データ(サンプル):MGDSST.MandA2012Leg3.track.txt
2012 06 16 14 50 26. 31.0200 126. 51.7510 300.632 299.933 300.726
2012 06 16 14 51 26. 30.8950 126. 51.5490 300.637 299.940 300.731
2012 06 16 14 52 26. 30.7710 126. 51.3470 300.641 299.946 300.735
2012 06 16 14 53 26. 30.6480 126. 51.1440 300.645 299.953 300.740
2012 06 16 14 54 26. 30.5230 126. 50.9430 300.649 299.960 300.744
2012 06 16 14 55 26. 30.3980 126. 50.7430 300.654 299.966 300.749
2012 06 16 14 56 26. 30.2750 126. 50.5440 300.658 299.973 300.754
2012 06 16 14 57 26. 30.1520 126. 50.3430 300.663 299.980 300.758
2012 06 16 14 58 26. 30.0290 126. 50.1440 300.667 299.987 300.763
2012 06 16 14 59 26. 29.9070 126. 49.9420 300.672 299.993 300.768