気象庁MGDSSTデータのバーチャル・サンプリング (航路上のデータを抽出する)

実行用シェルスクリプト: 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

&para

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