MSMデータとMGDSSTデータから潜熱・顕熱フラックスの平面図を作図する

Fortranプログラムのコンパイル

$ make clean; make

rm -rf ../obj/msm.heatflux.map.o ../obj/sh_lh_kagimoto.o ../obj/*.

real 0m0.006s

user 0m0.003s

sys 0m0.004s

if [ ! -d ../obj ]; then \

mkdir -p ../obj ; \

fi

ifort -c -traceback -module ../obj -assume byterecl -c -o ../olux.map.o msm.heatflux.map.f90

ifort -c -traceback -module ../obj -assume byterecl -c -o ../oimoto.o sh_lh_kagimoto.f90

ifort -o msm.heatflux.map ../obj/msm.heatflux.map.o ../obj/sh_lh_kadule ../obj -L/usr/local/lib -lnetcdf -assume byterecl

real 0m0.874s

user 0m0.730s

sys 0m0.144s

プログラムの実行

$ msm.heatflux.map.run.sh

Shell script, msm.heatflux.map.run.sh starts.

Input: /work4/DATA/MSM/MSM-S/2012/0614.nc

Input: /work2/kunoki/to_manda_sensei/Tools/MGDSST.on.MSM/output 0614.bin

Reading u

Done!

Reading v

Done!

Reading temp

Done!

Reading sp

Done!

Reading rh

Done!

Reading lon

Done!

Reading lat

Done!

Reading time

Done!

Reading SST file ...

Done.

Output: output/heatflux.msm.20120614.bin

Done ./msm.heatflux.map.run.sh

GrADSによる作図

ga-> open heatflux.msm.20120614.ctl

Scanning description file: heatflux.msm.20120614.ctl

Data file heatflux.msm.%y4%m2%d2.bin is open as file 1

LON set to 120 150

LAT set to 22.4 47.6

LEV set to 0 0

Time values set: 2012:6:14:0 2012:6:14:0

E set to 1 1

ga-> d lhf

Fortranプログラム実行用シェルスクリプト

#!/bin/bash

# Description:

#

# Author: kunoki

#

# Host: aofd165.fish.nagasaki-u.ac.jp

# Directory: /work2/kunoki/to_manda_sensei/Tools/MSM.HeatFlux.Map

#

# Revision history:

# This file is created by /usr/local/mybin/nbscr.sh at 13:35 on 10-24-2014.

echo "Shell script, $(basename $0) starts."

echo

exe=msm.heatflux.map

iyyyy=2012

imm=06

idd=14

yyyymmdd=${iyyyy}${imm}${idd}

indir=/work4/DATA/MSM/MSM-S/${iyyyy}

in=${imm}${idd}.nc

infle=${indir}/$in

if [ ! -f $infle ]; then

echo Error in $: No such file, $infle

exit 1

fi

indir_sst=/work2/kunoki/to_manda_sensei/Tools/MGDSST.on.MSM/output

in_sst=mgdsst.msm.${yyyymmdd}.bin

infle_sst=${indir_sst}/$in_sst

if [ ! -f $infle_sst ]; then

echo Error in $: No such file, $infle_sst

exit 1

fi

outdir=output

out=heatflux.msm.${yyyymmdd}.bin

mkdir -vp $outdir

ofle=${outdir}/${out}

im=481

jm=505

lm=24

undef=-1.0000000E+35

namelist=${exe}.namelist.txt

cat <<EOF > $namelist

&para

infle="${infle}",

im=${im}

jm=${jm}

lm=${lm}

infle_sst="${infle_sst}",

iyyyy=${iyyyy},

imm=${imm},

idd=${idd},

undef=${undef},

ofle="${ofle}"

&end

EOF

${exe} < $namelist

if [ ${imm} -eq 1 ]; then

mmm=Jan

elif [ ${imm} -eq 2 ]; then

mmm=Feb

elif [ ${imm} -eq 3 ]; then

mmm=Mar

elif [ ${imm} -eq 4 ]; then

mmm=Apr

elif [ ${imm} -eq 5 ]; then

mmm=May

elif [ ${imm} -eq 6 ]; then

mmm=Jun

elif [ ${imm} -eq 7 ]; then

mmm=Jul

elif [ ${imm} -eq 8 ]; then

mmm=Aug

elif [ ${imm} -eq 9 ]; then

mmm=Sep

elif [ ${imm} -eq 10 ]; then

mmm=Nov

elif [ ${imm} -eq 11 ]; then

mmm=Dec

fi

ctl=${outdir}/heatflux.msm.${yyyymmdd}.ctl

cat <<EOF >$ctl

dset ^heatflux.msm.%y4%m2%d2.bin

title heat flux

undef ${undef}

options template little_endian yrev

xdef ${im} linear 120.0 0.0625

ydef ${jm} linear 22.4 0.0500

zdef 1 linear 0.0 1

tdef 24 linear 00z${idd}${mmm}${iyyyy} 1hr

vars 8

lhf 1 0 Latent Heat Flux in W/m2

shf 1 0 Sensible Heat Flux in W/m2

u 1 0 10m-eastward wind m/s

v 1 0 10m-northward wind m/s

temp 1 0 surface air temp in deg. C

sst 1 0 sea surface temp in deg. C

rh 1 0 relative humidity in %

sp 1 0 surface pressure in hPa

endvars

EOF

echo "Done $0"

echo

Fortranプログラムコンパイル用makefile

#

#

#j Fortran77とFortran90/95の混成プログラムのコンパイル用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=-traceback #-CB -fpe0

#j コンパイルオプション

FFLAGS=-module ${OBJDIR} -assume byterecl

#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 -assume byterecl

#j ターゲット名(最終的に作りたい実行ファイル名)

TARGET1=msm.heatflux.map

TARGET2=

TARGETS=$(TARGET1) $(TARGET2)

#j MOD (Fortran90以降のモジュールファイル)

MOD=

#j SRC7 (Fortran77のソースファイル)

SRC7=

#j SRC9 (Fortran90のソースファイル)

SRC9=msm.heatflux.map.f90 \

sh_lh_kagimoto.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プログラム

メインルーチン: msm_heatflux_map.f90

program msm_heatflux_map

! Description:

!

! Author: kunoki

!

! Host: aofd165.fish.nagasaki-u.ac.jp

! Directory: /work2/kunoki/to_manda_sensei/Tools/MSM.HeatFlux.Map

!

! Revision history:

! Created by /usr/local/mybin/nff.sh at 11:24 on 10-24-2014.

! References

! http://www.rain.hyarc.nagoya-u.ac.jp/laboratory/OB/shimizu/HTML/netcdf.html

! Tips: NetCDF. http://iprc.soest.hawaii.edu/users/sasakiyo/Japanese/misc/tips_netcdf.html

! NetCDF Tips. http://www.sci.hokudai.ac.jp/grp/poc/top/software/other/netcdf_tips/index.htm

include '/usr/local/include/netcdf.inc'

! ncid:ファイルのID番号、 varid: 変 数のID番号

integer ncid,varid

character infle*500,infle_sst*500, ofle*500, ofle2*500, ofle2_base*500, ofle3*500

character ofle_tmp*500, ofle_tmp2*500

character var*50,day*50,largeday*50,largeday2*50

real,allocatable :: lon(:),lat(:),p(:),time(:)

integer(2),allocatable :: uin(:,:,:), vin(:,:,:),tempin(:,:,:), &

& spin(:,:,:)

integer(2),allocatable :: rhin(:,:,:)

integer,allocatable :: dimids(:),nshape(:)

character*70 err_message

character*100 dimname !dimnameは各次元の名前。

real,parameter::RMISS=-999.9 !欠損値

real,allocatable :: u(:,:,:),v(:,:,:),temp(:,:,:),sp(:,:,:),rh(:,:,:)

real,allocatable :: sst(:,:)

real,allocatable :: qh(:,:,:), qe(:,:,:)

! qh: Sensible Heat flux

! qe: Latent Heat flux

real :: undef

real::scale_factor, add_offset

character arg1*10,arg2*10,arg3*10

character month*3

integer stat,sta

integer wlength

namelist /para/infle,im,jm,lm,infle_sst,iyyyy,imm,idd,undef,ofle

read(*,nml=para)

allocate(lon(im),lat(jm),time(lm),u(im,jm,lm),v(im,jm,lm),&

& temp(im,jm,lm),sp(im,jm,lm))

allocate(uin(im,jm,lm),vin(im,jm,lm),tempin(im,jm,lm), &

& spin(im,jm,lm))

allocate(rhin(im,jm,lm),rh(im,jm,lm))

allocate(sst(im,jm))

allocate(qh(im,jm,lm),qe(im,jm,lm))

print '(A,A)','Input: ',trim(infle)

print '(A,A)','Input: ',trim(infle_sst)

print *

sta = nf_open(infle, nf_nowrite, ncid) ! ファイルのopenとNetCDF ID(ncid)の取得

if(sta.ne.0) then

err_message = nf_strerror(sta)

write(6,*) err_message

stop "Error : Cannot read input file."

endif

var='u'

print '(A,A)','Reading ',var(1:lnblnk(var))

sta = nf_inq_varid(ncid,var,varid)

sta = nf_get_att_real(ncid,varid,'scale_factor',scale_u)

sta = nf_get_att_real(ncid,varid,'add_offset',offset_u)

stat = nf_inq_varndims(ncid, varid, ndims)

allocate(dimids(ndims))

stat = nf_inq_vardimid(ncid, varid, dimids)

allocate(nshape(ndims))

do i=1,ndims

stat = nf_inq_dimlen(ncid, dimids(i), nshape(i))

enddo !i

do i=1,ndims

stat = nf_inq_dimname(ncid, dimids(i), dimname)

enddo

stat = nf_get_var_INT2(ncid, varid, uin)

deallocate(dimids,nshape)

print '(A)','Done!'

print *

var='v'

print '(A,A)','Reading ',var(1:lnblnk(var))

sta = nf_inq_varid(ncid,var,varid)

sta = nf_get_att_real(ncid,varid,'scale_factor',scale_v)

sta = nf_get_att_real(ncid,varid,'add_offset',offset_v)

stat = nf_inq_varndims(ncid, varid, ndims)

allocate(dimids(ndims))

stat = nf_inq_vardimid(ncid, varid, dimids)

allocate(nshape(ndims))

do i=1,ndims

stat = nf_inq_dimlen(ncid, dimids(i), nshape(i))

enddo !i

do i=1,ndims

stat = nf_inq_dimname(ncid, dimids(i), dimname)

enddo

stat = nf_get_var_INT2(ncid, varid, vin)

deallocate(dimids,nshape)

print '(A)','Done!'

print *

var='temp'

print '(A,A)','Reading ',var(1:lnblnk(var))

sta = nf_inq_varid(ncid,var,varid)

sta = nf_get_att_real(ncid,varid,'scale_factor',scale_temp)

sta = nf_get_att_real(ncid,varid,'add_offset',offset_temp)

stat = nf_inq_varndims(ncid, varid, ndims)

allocate(dimids(ndims))

stat = nf_inq_vardimid(ncid, varid, dimids)

allocate(nshape(ndims))

do i=1,ndims

stat = nf_inq_dimlen(ncid, dimids(i), nshape(i))

enddo !i

do i=1,ndims

stat = nf_inq_dimname(ncid, dimids(i), dimname)

enddo

stat = nf_get_var_INT2(ncid, varid, tempin)

deallocate(dimids,nshape)

print '(A)','Done!'

print *

var='sp'

print '(A,A)','Reading ',var(1:lnblnk(var))

sta = nf_inq_varid(ncid,var,varid)

sta = nf_get_att_real(ncid,varid,'scale_factor',scale_sp)

sta = nf_get_att_real(ncid,varid,'add_offset',offset_sp)

stat = nf_inq_varndims(ncid, varid, ndims)

allocate(dimids(ndims))

stat = nf_inq_vardimid(ncid, varid, dimids)

allocate(nshape(ndims))

do i=1,ndims

stat = nf_inq_dimlen(ncid, dimids(i), nshape(i))

enddo !i

do i=1,ndims

stat = nf_inq_dimname(ncid, dimids(i), dimname)

enddo

stat = nf_get_var_INT2(ncid, varid, spin)

deallocate(dimids,nshape)

print '(A)','Done!'

print *

var='rh'

print '(A,A)','Reading ',var(1:lnblnk(var))

sta = nf_inq_varid(ncid,var,varid)

sta = nf_get_att_real(ncid,varid,'scale_factor',scale_rh)

sta = nf_get_att_real(ncid,varid,'add_offset',offset_rh)

stat = nf_inq_varndims(ncid, varid, ndims)

allocate(dimids(ndims))

stat = nf_inq_vardimid(ncid, varid, dimids)

allocate(nshape(ndims))

do i=1,ndims

stat = nf_inq_dimlen(ncid, dimids(i), nshape(i))

enddo !i

do i=1,ndims

stat = nf_inq_dimname(ncid, dimids(i), dimname)

enddo

stat = nf_get_var_INT2(ncid, varid, rhin)

deallocate(dimids,nshape)

print '(A)','Done!'

print *

var='lon'

print '(A,A)','Reading ',var(1:lnblnk(var))

sta = nf_inq_varid(ncid,var,varid)

stat = nf_get_var_real(ncid, varid, lon)

print '(A)','Done!'

print *

var='lat'

print '(A,A)','Reading ',var(1:lnblnk(var))

sta = nf_inq_varid(ncid,var,varid)

stat = nf_get_var_real(ncid, varid, lat)

print '(A)','Done!'

print *

var='time'

print '(A,A)','Reading ',var(1:lnblnk(var))

sta = nf_inq_varid(ncid,var,varid)

stat = nf_get_var_real(ncid, varid, time)

print '(A)','Done!'

print *

print '(A)','Reading SST file ...'

inquire(iolength=wlength) sst(1,1)

siz = jm * im * wlength

open(11,file=infle_sst,access='direct',form='unformatted', &

& & action="read",recl=siz)

read(11,rec=1) ((sst(i,j),i=1,im),j=1,jm)

close(11)

print '(A)','Done.'

print *

do i=1,im

do j=1,jm

do l=1,lm

u(i,j,l)=uin(i,j,l)*scale_u + offset_u

v(i,j,l)=vin(i,j,l)*scale_v + offset_v

temp(i,j,l)=tempin(i,j,l)*scale_temp + offset_temp - 273.15

sp(i,j,l)=(spin(i,j,l)*scale_sp + offset_sp)/100.0

rh(i,j,l)=rhin(i,j,l)*scale_rh + offset_rh

if(u(i,j,l)>=-100.0 .and. v(i,j,l)>=-100.0 .and. temp(i,j,l)>=-50.0 &

.and. sp(i,j,l)>=800.0 .and. rh(i,j,l)>=0.0 .and. sst(i,j)>=-100)then

v10=sqrt(u(i,j,l)**2 + v(i,j,l)**2)

call sh_lh_kagimoto(sp(i,j,l),sst(i,j),temp(i,j,l),rh(i,j,l),v10,&

& qh(i,j,l),qe(i,j,l))

else

qh(i,j,l)=undef

qe(i,j,l)=undef

endif

enddo !l

enddo !i

enddo !j

print '(A,A)','Output: ',trim(ofle)

open(21,file=ofle,form="unformatted",access='direct',&

& recl=im*jm*4)

irec=1

do l=1,lm

write(21,rec=irec) ((qe(i,j,l),i=1,im),j=1,jm)

irec=irec+1

write(21,rec=irec) ((qh(i,j,l),i=1,im),j=1,jm)

irec=irec+1

write(21,rec=irec) ((u(i,j,l),i=1,im),j=1,jm)

irec=irec+1

write(21,rec=irec) ((v(i,j,l),i=1,im),j=1,jm)

irec=irec+1

write(21,rec=irec) ((temp(i,j,l),i=1,im),j=1,jm)

irec=irec+1

write(21,rec=irec) ((sst(i,j),i=1,im),j=1,jm)

irec=irec+1

write(21,rec=irec) ((rh(i,j,l),i=1,im),j=1,jm)

irec=irec+1

write(21,rec=irec) ((sp(i,j,l),i=1,im),j=1,jm)

irec=irec+1

enddo !L

close(21)

end program msm_heatflux_map

サブルーチン: sh_lh_kagimoto.f90

subroutine sh_lh_kagimoto(slp,sst,t,rh,v10,qh,qe)

! Description:

!

! Author: am

!

! Host: aofd30

! Directory: /work2/am/12.Work11/23.Recomp_HeatFlux/23.4Fluxes_Kagimoto/src

!

! Revision history:

! This file is created by /usr/local/mybin/nff.sh at 08:50 on 11-14-2011.

! use

! implicit none

real,intent(in):: slp,sst,t,rh,v10

! t: SAT in deg. C

! sst: SST in deg. C

! slp: Sea Level Pressure in hPa

! rh : Relative Humidity in %

! v10: Wind speed at a height of 10 m in m/s

real,intent(out)::qh,qe

real ae(5)/0., 0.969, 1.18, 1.196, 1.68/

real ah(5)/0., 0.927, 1.15, 1.17, 1.652/

real be(5)/1.23, 0.0521, 0.01, 0.008, -0.016/

real bh(5)/1.185, 0.0546, 0.01, 0.0075, -0.017/

real ce(5)/0., 0, 0, -0.0004, 0./

real ch(5)/0., 0, 0, -0.00045, 0./

real pe(5)/-0.16, 1., 1., 1., 1./

real ph(5)/-0.157, 1., 1., 1., 1./

real,parameter :: dens=1.225 !kg/m3

! z0: roughness depth

real,parameter :: z0=1. ![m]

! z0=1: typical value as used in Skyllingstad et al. (JPO, 1999)

! (p. 730L, Noh et al., JPO, 2004)

real es,essea,q,qs !,dens

real v ! Wind speed at 2m above the surface

real chn,cln,chg,clg,s0,s

real LH,cp

parameter(LH=2501.,cp=1005.)

! write(*,'(a)')'Subroutine: sh_lh_kagimoto'

! write(*,*)''

z1=10.0 ![m] height of the wind measurement

z=2.

! v=v10*log(z/z0)/log(z1/z0)

v=v10

! print *,'v=' ,v

if (v.le.2.2) l=1

if (v.gt.2.2 .and. v.le.5.0) l=2

if (v.gt.5.0 .and. v.le.8.0) l=3

if (v.gt.8.0 .and. v.le.25.0) l=4

if (v.gt.25.0 .and. v.le.50.0) l=5

chn=0.001*(ah(l) + bh(l)*v**ph(l) + ch(l)*(v-8.)**2)

cln=0.001*(ae(l) + be(l)*v**pe(l) + ce(l)*(v-8.)**2)

s0=(sst-t)/(v**2)

s=s0*abs(s0)/(abs(s0)+0.01)

if(s.lt. -3.3) then

chg=0.0

clg=0.0

else if(s .lt. 0.0) then

chg=chn*(0.1+0.03*s+0.9*exp(4.8*s))

clg=cln*(0.1+0.03*s+0.9*exp(4.8*s))

else

chg=chn*(1.0+0.63*sqrt(s))

clg=cln*(1.0+0.63*sqrt(s))

endif

es=6.11*exp((LH/0.4615)*(1.0/273.2-1.0/(t+273.15)))

essea=6.11*exp((LH/0.4615)*(1.0/273.2-1.0/(sst+273.15)))*100/rh

q=0.622*es/slp

qs=0.622*essea/slp

! dens=slp/(2.87*(1+0.61*q)*(t+273.15))

! print *,'dens=' ,dens

qh=dens*cp*chg*v*(sst-t)

qe=dens*LH*1000*clg*v*(qs-q)

end subroutine sh_lh_kagimoto