[2017年 1月 12日 木曜日 14:15:59 JST]
[~/2016.PolarLow/Ship.Surface.Flux]
[am@aofd165]
$ srcdump.sh flux.PL1101.sh flux.PL1101.f90 flux.PL1101.plot.gmt.sh flux.PL1101.plot.bowen.gmt.sh gmtpar.sh note.sh>tmp.txt
------------------------------
List of the following files:
------------------------------
flux.PL1101.sh
flux.PL1101.f90
flux.PL1101.plot.gmt.sh
flux.PL1101.plot.bowen.gmt.sh
gmtpar.sh
note.sh
------------------------------
Machine info
------------------------------
aofd165.bio.mie-u.ac.jp
/work1/am/2016.PolarLow/Ship.Surface.Flux
Thu Jan 12 14:17:08 JST 2017
======================
flux.PL1101.sh
======================
#!/bin/bash
nml=$(basename $0 .sh).namelist.txt
src=$(basename $0 .sh).f90
exe=$(basename $0 .sh).exe
err=$(basename $0 .sh).err
fort=ifort
opt="-CB -traceback -fpe0 -assume byterecl"
temp=$(basename $src .f90)
prog_name=$(echo $temp| sed -e 's/\./_/g')
echo Compiling ${src} ...
${fort} ${opt} ${src} -o ${exe}
if [ $? -ne 0 ]; then
echo
echo Compile error!
echo
echo Terminated.
echo
exit 1
fi
echo "Done Compile ${src}."
ls -lh ${exe}
echo
input="PL1101.meteo.for.flux.txt"
nw_half=2
input_sst="PL1101.SST.for.flux.txt"
nline=$(wc $input_sst | awk '{print $1}')
nd_sst=$(expr $nline / 10)
rmiss=-999.0
output="PL1101.FLUX.txt"
output_check="PL1101.FLUX.CHECK.txt"
cat <<EOF>$nml
¶
infle="${input}",
infle_sst="${input_sst}",
ofle="${output}",
ofle_check="${output_check}",
nw_half=${nw_half}
rmiss=${rmiss},
nd_sst=${nd_sst}
&end
EOF
${exe} <${nml}
exit
----------------------
End of flux.PL1101.sh
----------------------
======================
flux.PL1101.f90
======================
program flux_PL1101
character(len=500 )::infle,infle_sst, ofle, ofle_check
character(len=2000)::strm
character(len=2000),allocatable::comment(:)
character(len=19),allocatable ::datetime(:)
character(len=19)::sdatetime
character(len=8)::sdate
character(len=6)::stime
real::elptime
integer,parameter::nvar=10
real,parameter::r2d=180.0/3.14159265358979
real,allocatable:: var_in(:,:), var_out(:,:)
real,allocatable:: v1_in(:), v1_out(:)
real,allocatable:: ws_out(:), wd_out(:)
real(8),allocatable:: elptime_met(:)
character(len=19),allocatable ::datetime_sst(:)
real(8),allocatable:: elptime_sst(:)
real, allocatable:: sst(:), depth(:)
real,allocatable:: lond(:), lonm(:), latd(:), latm(:)
integer,allocatable::vld(:)
character(len=300),allocatable:: infle_t(:)
integer,allocatable::yyyy(:),mm(:),dd(:),hh(:),nn(:),ss(:)
integer ::syyyy, smm, sdd, shh, snn, sss
real(8),allocatable::elptime_nearest(:)
integer,allocatable::ndx(:)
real(8)::dtime_min, d_elptime
real,allocatable::qe(:),qh(:),ev(:)
real rmiss
namelist /para/infle,infle_sst,ofle,ofle_check,nw_half,rmiss,nd_sst
print *
print '(A)','Read nemelist file.'
read(*,nml=para)
print '(A)','Done.'
print *
print '(A)','Read meteorological data.'
open(11,file=infle,action="read")
n=0
nc=0
count_vaild_data: do
read(11,'(A)',iostat=ios)strm
if(ios<0)exit
if(strm(1:1) == "#")then
nc=nc+1
cycle
else
n=n+1
endif
enddo count_vaild_data
allocate(datetime(n),elptime_met(n))
allocate(var_in(nvar,n), var_out(nvar,n))
allocate(v1_in(n), v1_out(n))
allocate(ws_out(n), wd_out(n))
allocate(yyyy(n),mm(n),dd(n),hh(n),nn(n),ss(n))
nt=n
allocate(comment(nc))
rewind(11)
n=0
nc=0
skip_comment: do
read(11,'(A)',iostat=ios)strm
if(ios<0)exit
if(strm(1:1) == "#")then
nc=nc+1
comment(nc)=strm
cycle
else
n=n+1
read(strm(1:19),'(A19)') datetime(n)
read(strm(20:) ,* )(var_in(i,n),i=1,nvar)
end if
enddo skip_comment
ntcheck=n
print '(A,2i5)',"Number of input data lines = ",nt,ntcheck
close(11)
print '(A)','Done.'
print *
print '(A)','Read SST data.'
allocate(infle_t(nd_sst))
allocate(datetime_sst(nd_sst))
allocate(elptime_sst(nd_sst),sst(nd_sst))
allocate(vld(nd_sst))
allocate(lond(nd_sst),lonm(nd_sst),latd(nd_sst),latm(nd_sst))
allocate(depth(nd_sst))
allocate(elptime_nearest(nd_sst), ndx(nd_sst))
open(12,file=infle_sst,action="read")
do n=1,nd_sst
read(12,'(A)') strm
read(strm(3:),'(A)') infle_t(n)
read(12,'(A)') strm
read(strm(10:),*) lond(n),lonm(n), latd(n), latm(n)
read(12,'(A)') strm
read(strm(9:),*) depth(n)
read(12,'(A)') strm
read(strm(14:),'(A)') sdatetime
read(12,'(A)') strm
read(strm(13:),'(A)') datetime_sst(n)
read(12,'(A)') strm
read(strm(7:),*) dz
read(12,'(A)') strm
read(strm(7:),*) vld(n)
read(12,'(A)') strm
read(strm(8:),*) zbtm
read(12,'(A)') strm
!
read(12,'(A)') strm
read(strm,*) elptime_sst(n), sst(n)
! print *, elptime_sst(n)
end do !n
print '(A)','Done.'
print *
print '(A)','Set origin for computing the elapsed time'
read(sdatetime(1:4), '(i4.4)')syyyy
read(sdatetime(6:7), '(i2.2)')smm
read(sdatetime(9:10),'(i2.2)')sdd
read(sdatetime(12:13),'(i2.2)')shh
read(sdatetime(15:16),'(i2.2)')snn
read(sdatetime(18:19),'(i2.2)')sss
print *,syyyy,smm, sdd,shh,snn,sss
print '(A)','Done.'
print *
print '(A)','Compute elapsed time'
do n=1,nt
read(datetime(n)(1:4), '(i4.4)')yyyy(n)
read(datetime(n)(6:7), '(i2.2)')mm(n)
read(datetime(n)(9:10),'(i2.2)')dd(n)
read(datetime(n)(12:13),'(i2.2)')hh(n)
read(datetime(n)(15:16),'(i2.2)')nn(n)
read(datetime(n)(18:19),'(i2.2)')ss(n)
call eltime(syyyy, smm, sdd, shh, snn, sss, &
yyyy(n), mm(n), dd(n), hh(n), nn(n), ss(n), &
elptime_met(n))
!print *,syyyy, smm, sdd, shh, smm, sss
!print *,yyyy(n), mm(n), dd(n), hh(n), mm(n), ss(n)
!print *,elptime_met(n)
enddo !n
print '(A)','Done.'
print *
print '(A)','Find the nearest date/time for flux computation'
do k=1,nd_sst
dtime_min=9.D9
ndx(k)=1
do n=1,nt
d_elptime = abs(elptime_sst(k)-elptime_met(n))
if(d_elptime < dtime_min)then
ndx(k)=n
dtime_min = d_elptime
end if
end do !n
elptime_nearest(k)=elptime_met(ndx(k))
! print *, datetime(ndx(k))," ", datetime_sst(k)
! print *, elptime_met(ndx(k))," ", elptime_sst(k)
end do !k
print '(A)','Done.'
print *
print '(A)','Low-pass filtering for met data'
do i=1,nvar
do n=1,nt
v1_in(n)=var_in(i,n)
v1_out(n)=rmiss
end do !n
call boxcar(nt, v1_in, nw_half, v1_out, rmiss)
do n=1,nt
var_out(i,n)=v1_out(n)
end do !n
end do !i
print '(A)','Done.'
print *
print '(A)','Compute wind direction of low-passed data'
do n=1,nt
if(var_out(8,n)/=rmiss .and. var_out(9,n)/=rmiss)then
ws_out(n)=sqrt(var_out(8,n)**2 + var_out(9,n)**2)
wd_out(n)=atan2(var_out(9,n),var_out(8,n))*r2d
else
ws_out(n)=rmiss
wd_out(n)=rmiss
endif
enddo !n
var_out(6,:) =ws_out(:)
var_out(7,:)=wd_out(:)
print '(A)','Done.'
print *
print '(A)','Compute sensible and latent heat flux'
allocate(qh(nd_sst), qe(nd_sst), ev(nd_sst))
! prs [Pa]
! sstk [K]
! sat [K]
! rh [%]
! u10 [m/s]
! v10 [m/s]
do k=1,nd_sst
prs= var_out(3,ndx(k))*100.0
sat= var_out(1,ndx(k))+273.15
sstk=sst(k) +273.15
u10= var_out(8,ndx(k))
v10= var_out(9,ndx(k))
rh = var_out(2,ndx(k))
call flux_kondo(prs, sat, rh, u10, v10, sstk, rmiss, qh(k), qe(k), &
& ev(k))
end do !n
print '(A)','Done.'
print *
print '(A)','Output check data'
open(21,file=ofle_check)
write(21,'(A,A)')'# Input: ',trim(infle)
write(21,'(A)')'# Date/Time SAT RH SLP Lat&
& Lon ws wd we wn q sst'
write(21,'(A)')'# degC % hPa &
° N deg E m/s deg m/s m/s g/kg degC'
do k=1,nd_sst
write(21,'(A19,2f8.2,f9.2,2f10.4,2f10.2,2f9.3,f9.2, f10.2)')&
datetime_sst(k),(var_out(i,ndx(k)),i=1,nvar), sst(k)
enddo !n
print '(A)','Done.'
print *
print '(A)','Output flux data'
open(21,file=ofle)
write(21,'(A,A)')'# Input: ',trim(infle)
write(21,'(A)')'# Positive upward'
write(21,'(A)')'# Date/Time SHF LHF E'
write(21,'(A)')'# W/m2 W/m2 mm/h'
do k=1,nd_sst
write(21,'(A19,2f9.2,f10.4)')&
datetime_sst(k),qh(k),qe(k),ev(k)
enddo !n
close(21)
print '(A)','Done.'
print *
print '(A,A)','Input: ',trim(infle)
print '(A,A)','Input: ',trim(infle_sst)
print '(A,A)','Output: ',trim(ofle)
print '(A,A)','Output: ',trim(ofle_check)
print *
end program flux_PL1101
subroutine flux_kondo(prs, sat, rh, u10, v10, sstk, rmiss, qh, qe, ev)
!
! Description: Surface heat flux by Kondo
!
! Author: am
!
! Host: aofd165.fish.nagasaki-u.ac.jp
! Directory: /work1/am/Sato.Kodama.Manda2015/MSM.MGDSST.FLUX.Trajectory
!
real,intent(in):: prs, sat, rh, u10, v10, sstk, rmiss
real,intent(out)::qh, qe, ev
! slp [Pa]
! sstk [K]
! sat [K]
! rh [%]
! u10 [m/s]
! v10 [m/s]
! q specific humidity [kg/kg]
! qh SHF [W/m2]
! qe: LHF [W/m2]
! ev: Evapolation [mm/h]
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 es,essea,qs,dens
real chn,cln,chg,clg,s0,s
parameter(LH=2500,cp=1004)
if(prs==rmiss .or. sat==rmiss .or. rh==rmiss .or. u10==rmiss .or. &
v10==rmiss .or. sstk==rmiss)then
qh=rmiss
qe=rmiss
ev=rmiss
return
endif
slp=prs/100.0 !Pa -> hPa
sst=sstk-273.15 !K -> degC
t=sat-273.15
v=sqrt(u10**2+v10**2)
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))
! write(6,*) 'dens,s,chg,clg,es,essea,q,qs'
! write(6,*) dens,q,qs
qh=dens*cp*chg*v*(sst-t)
qe=dens*LH*1000*clg*v*(qs-q)
ev=qe/(LH*1000.0)*3600.0 !mm/h
end subroutine flux_kondo
subroutine eltime(yyyy1,mo1,dd1,hh1,mi1,ss1, &
yyyy2,mo2,dd2,hh2,mi2,ss2, &
elapsed_time )
integer,intent(in)::yyyy1,mo1,dd1,hh1,mi1,ss1
integer,intent(in)::yyyy2,mo2,dd2,hh2,mi2,ss2
real(8),intent(inout)::elapsed_time
integer::julian1, julian2
real(8)::timeh1, timeh2
real(8),parameter::d2h=24.d0, m2h=1.d0/60.d0, s2h=1.d0/3600.d0
julian1=julday(yyyy1,mo1,dd1)
julian2=julday(yyyy2,mo2,dd2)
timeh1=dble(julian1)*d2h + dble(hh1) + dble(mi1)*m2h + dble(ss1)*s2h
timeh2=dble(julian2)*d2h + dble(hh2) + dble(mi2)*m2h + dble(ss2)*s2h
elapsed_time=timeh2-timeh1
! print *,timeh2, timeh1
return
end subroutine eltime
subroutine caldat(JULIAN,IYYY,MONTH,DD)
INTEGER,intent(in):: JULIAN
integer,intent(out)::IYYY,MONTH,DD
! NAME IN/OUT DESCRIPTION
!
! JULIAN I THE JULIAN DAY
! IYYY O THE YEAR
! MONTH O THE MONTH (1 TO 12)
! DD O THE DAY OF THE MONTH
!
INTEGER IGREG
PARAMETER (IGREG=2299161)
INTEGER JALPHA,JA,JB,JC,JD,JE
!
IF (JULIAN .GE. IGREG) THEN
JALPHA = INT(((JULIAN - 1867216) - 0.25)/36524.25)
JA = JULIAN + 1 + JALPHA - INT(0.25*JALPHA)
ELSE
JA = JULIAN
ENDIF
JB = JA + 1524
JC = INT(6680. + ((JB - 2439870) - 122.1)/365.25)
JD = 365*JC + INT(0.25*JC)
JE = INT((JB - JD)/30.6001)
DD = JB - JD - INT(30.6001*JE)
MONTH = JE - 1
IF (MONTH .GT. 12) MONTH = MONTH - 12
IYYY = JC - 4715
IF (MONTH .GT. 2) IYYY = IYYY - 1
IF (IYYY .LE. 0) IYYY = IYYY - 1
RETURN
end subroutine caldat
INTEGER FUNCTION JULDAY(IYYY,MONTH,DD)
INTEGER:: IYYY,MONTH,DD
!
! NAME IN/OUT DESCRIPTION
!
! IYYY I YEAR
! MONTH I MONTH (1 TO 12)
! DD I DAY OF MONTH
! JULDAY O JULIAN DAY
!
INTEGER IGREG
PARAMETER (IGREG = 15 + 31*(10 + 12*1582))
INTEGER JY,JM,JA
!
IF (IYYY .LT. 0) IYYY = IYYY + 1
IF (MONTH .GT. 2) THEN
JY = IYYY
JM = MONTH + 1
ELSE
JY = IYYY - 1
JM = MONTH + 13
ENDIF
JULDAY = INT(365.25*JY) + INT(30.6001*JM) + DD + 1720995
IF (DD + 31*(MONTH + 12*IYYY) .GE. IGREG) THEN
JA = INT(0.01*JY)
JULDAY = JULDAY + 2 - JA + INT(0.25*JA)
ENDIF
RETURN
end function julday
subroutine usage()
print '(A)','Usage: eltime yyyymmdd1 hhmmss1 yyyymodd2 hhmiss2'
print '(A)','yyyymmdd: 8-digit number representing year, month, day'
print '(A)','hhmmss: 6-digit number representing hour, minute, and second.'
end subroutine
subroutine boxcar(nt, v_in, nw_half, v_out, rmiss)
integer,intent(in)::nt, nw_half
real,intent(in):: v_in(nt)
real,intent(inout)::v_out(nt)
real,intent(in)::rmiss
integer,parameter::zero_i=0
real, parameter::zero_r=0.0
do n = 1 + nw_half, nt - hw_half -2
n_valid=zero_i
sum=zero_r
do i = -nw_half, nw_half
if(v_in(n+i)/=rmiss)then
n_valid=n_valid+1
sum = sum + v_in(n+i)
end if
end do !i
if(n_valid > nw_half)then
v_out(n)=sum/float(n_valid)
else
v_out(n)=rmiss
endif
end do !n
return
end subroutine boxcar
----------------------
End of flux.PL1101.f90
----------------------
======================
flux.PL1101.plot.gmt.sh
======================
#!/bin/bash
# Description:
#
# Author: am
#
# Host: aofd165.bio.mie-u.ac.jp
# Directory: /work1/am/WRF.Post/NCL/PL1101/Fink.01.Dp
#
# Revision history:
# This file is created by /usr/local/mybin/ngmt.sh at 14:08 on 11-30-2016.
. ./gmtpar.sh
echo "Bash script $0 starts."
gmtset PLOT_CLOCK_FORMAT hh
gmtset PLOT_DATE_FORMAT "o dd"
gmtset TIME_FORMAT_SECONDARY abbreviated
gmtset ANOT_FONT_SIZE 12
gmtset LABEL_FONT_SIZE 14
gmtset ANNOT_FONT_SIZE_PRIMARY 10
gmtset ANNOT_FONT_SIZE_SECONDARY 12
#gmtset INPUT_DATE_FORMAT yyyy-mm-dd
#gmtset INPUT_CLOCK_FORMAT hh:mm:ss
range="0/50/68/80"
size="s25/90/3.5i/68"
anot=a10f5g5/5g5WSne
xrange=2011-01-19T00:00:00/2011-01-26T00:00:00
range1=${xrange}/-100/600
range2=${xrange}/975/1030
range3=${xrange}/-20/10
range4=${xrange}/50/100
range5=${xrange}/0/6
range6=${xrange}/0/20
size2="4T/1"
in0=PL1101.meteo.for.flux.lowpass.txt
in1=PL1101.FLUX.txt
in2=PL1101.FLUX.CHECK.txt
inlist="$in0 $in1 $in2"
for infile in $inlist; do
if [ ! -f $infile ]; then
echo Error in $0 : No such file, $infile
exit 1
fi
done
#figdir="./fig/"
#if [ ! -d ${figdir} ];then
# mkdir -p $figdir
#fi
#out=${figdir}$(basename $in .asc).ps
out=${figdir}$(basename $0 .sh)_$(basename $in1 .txt).ps
echo Input : $in0
echo Input : $in1
echo Input : $in2
echo Output : $out
pscoast -R$range -J$size -G200 -W1 -Di \
-X1.5 -Y7.5 -P -K > $out
awk '{if ($1 != "#" && $5 >-900.0 && $6 > -900.0) print $6, $5}' $in0 |\
psxy -R$range -J$size -Sc0.01 -G0 \
-O -K >>$out
awk '{if ($1 != "#" && substr($1,12,8) == "00:00:00") print $6,$5}' $in0 |\
psxy -R$range -J$size -B$anot -Sx0.2 -W5 \
-O -K>>$out
awk '{if ($1 != "#" && $2 != -999.00) print $1, $2}' $in1|\
psxy -R$range1 -JX$size2 \
-Bpa12hf6h/a200f200:"[Wm@+-2@+]":Wsne -Bsa1D/f200g100000 \
-Sc0.05 -W3 \
-Y-1.5 -O -K >> $out
awk '{if ($1 != "#" && $3 != -999.00) print $1, $3}' $in1|\
psxy -R$range1 -JX$size2 \
-St0.07 -W3 \
-O -K >> $out
awk '{if ($1 != "#" && $2 > -900) print $1, $2}' $in2|\
psxy -R$range3 -JX$size2 \
-Bpa12hf6h/a5f5:"T${sp}[${deg}C]":Wsne -Bsa1D/f5 \
-W5 \
-Y-1.3 -O -K >> $out
awk '{if ($1 != "#" && $12 > -900) print $1, $12}' $in2|\
psxy -R$range3 -JX$size2 \
-Sc0.05 \
-O -K >> $out
awk '{if ($1 != "#" && $3 > -900) print $1, $3}' $in2|\
psxy -R$range4 -JX$size2 -Bpa12hf6h/a20f10:"RH${sp}[%]":Wsne -Bsa1D/f10 -W5 \
-Y-1.3 -O -K >> $out
awk '{if ($1 != "#" && $11 > -900) print $1, $11}' $in2|\
psxy -R$range5 -JX$size2 -Bpa12hf6h/a1f1:"q${sp}[gkg@+-1@+]":Wsne -Bsa1D/f1 -W5 \
-Y-1.3 -O -K >> $out
awk '{if ($1 != "#" && $7 > -900) print $1, $7}' $in2|\
psxy -R$range6 -JX$size2 -Bpa12hf6h/a10f5:"W@-S@-${sp}[ms@+-1@+]":WSne -Bsa1D/f5 -W5 \
-Y-1.3 -O -K >> $out
in="$in1 $in2"
xoffset=
yoffset=8.8
comment=$in0
. ./note.sh
echo "Done $0"
----------------------
End of flux.PL1101.plot.gmt.sh
----------------------
======================
flux.PL1101.plot.bowen.gmt.sh
======================
#!/bin/bash
# Description:
#
# Author: am
#
# Host: aofd165.bio.mie-u.ac.jp
# Directory: /work1/am/WRF.Post/NCL/PL1101/Fink.01.Dp
#
# Revision history:
# This file is created by /usr/local/mybin/ngmt.sh at 14:08 on 11-30-2016.
. ./gmtpar.sh
echo "Bash script $0 starts."
gmtset PLOT_CLOCK_FORMAT hh
gmtset PLOT_DATE_FORMAT "o dd"
gmtset TIME_FORMAT_SECONDARY abbreviated
gmtset ANOT_FONT_SIZE 12
gmtset LABEL_FONT_SIZE 14
gmtset ANNOT_FONT_SIZE_PRIMARY 10
gmtset ANNOT_FONT_SIZE_SECONDARY 12
#gmtset INPUT_DATE_FORMAT yyyy-mm-dd
#gmtset INPUT_CLOCK_FORMAT hh:mm:ss
range="0/50/68/80"
size="s25/90/3.5i/68"
anot=a10f5g5/5g5WSne
xrange=2011-01-19T00:00:00/2011-01-26T00:00:00
range1=${xrange}/-100/1000
range2=${xrange}/-0.2/2
range2b=${xrange}/0.01/10.0
size2="4T/1"
size2b="4T/1l"
in0=PL1101.meteo.for.flux.lowpass.txt
in1=PL1101.FLUX.txt
in2=PL1101.FLUX.CHECK.txt
inlist="$in0 $in1 $in2"
for infile in $inlist; do
if [ ! -f $infile ]; then
echo Error in $0 : No such file, $infile
exit 1
fi
done
#figdir="./fig/"
#if [ ! -d ${figdir} ];then
# mkdir -p $figdir
#fi
#out=${figdir}$(basename $in .asc).ps
out=${figdir}$(basename $0 .sh)_$(basename $in1 .txt).BOWEN.ps
echo Input : $in0
echo Input : $in1
echo Output : $out
pscoast -R$range -J$size -G200 -W1 -Di \
-X1.5 -Y7.5 -P -K > $out
awk '{if ($1 != "#" && $5 >-900.0 && $6 > -900.0) print $6, $5}' $in0 |\
psxy -R$range -J$size -Sc0.01 -G0 \
-O -K >>$out
awk '{if ($1 != "#" && substr($1,12,8) == "00:00:00") print $6,$5}' $in0 |\
psxy -R$range -J$size -B$anot -Sx0.2 -W5 \
-O -K>>$out
awk '{if ($1 != "#" && $2 != -999.00) print $1, $2+$3}' $in1|\
psxy -R$range1 -JX$size2 \
-Bpa12hf6h/a500f250:"[Wm@+-2@+]":WSn -Bsa1D/f500g100000 \
-W5/${red} \
-Y-1.5 -O -K >> $out
awk '{if ($1 !="#" && $2 !=-999.00 && $3 !=-999.00) \
print $1, $2/$3}' $in1|\
psxy -R$range2 -JX${size2} \
-Bpa12hf6h/a0.5f.5:"B":E -Bsa1D/f200g1 \
-W3/${blue} \
-O -K >> $out
awk '{if ($1 != "#" && $2 != -999.00) print $1, $2+$3}' $in1|\
psxy -R$range1 -JX$size2 \
-Bpa12hf6h/a500f250:"[Wm@+-2@+]":WSn -Bsa1D/f500g100000 \
-W5/${red} \
-Y-2 -O -K >> $out
awk '{if ($1 !="#" && $2 !=-999.00 && $3 !=-999.00) \
print $1, $2/$3}' $in1|\
psxy -R$range2b -JX${size2b} \
-Bpa12hf6h/a0.5f.5:"B":E -Bsa1D/f200g100000 \
-W3/${blue} \
-O -K >> $out
in="$in0 $in1"
xoffset=
yoffset=5.5
comment=$in0
. ./note.sh
echo "Done $0"
----------------------
End of flux.PL1101.plot.bowen.gmt.sh
----------------------
======================
gmtpar.sh
======================
#
# GMTのパラメータの設定 (version 4以上対象）
#
# 長さの単位はインチとする
gmtset MEASURE_UNIT INCH
gmtset LABEL_FONT_SIZE 18
gmtset HEADER_FONT_SIZE 20
gmtset ANOT_FONT_SIZE 18
gmtset TICK_PEN 4
# 地図の縦横軸に縞々を使わない
gmtset BASEMAP_TYPE PLAIN
# 紙のサイズはA4
gmtset PAPER_MEDIA A4
# 地図の°の記号
gmtset DEGREE_SYMBOL = degree # Available for ver. 4 or later
# 空白文字の設定
sp="\040" # White space
# =の記号
eq="\075" # equal
# 温度の°の記号
deg="\260" #deg="\312" # degree symbol
# 色のRGB値を設定 (線や点に色をつけるときRGB値を直接指定するより分かりやすい)
red="255/0/0"
orange="255/165/0"
pink="255/192/203"
green="0/255/0"
blue="0/0/255"
midnightblue="25/25/112"
yellow="255/255/0"
gold="255/215/0"
purple="160/32/240"
magenta="255/0/255"
white="255/255/255"
# 線種を指定
dash="t15_5:15"
dot="t5_5:5"
dotdash="t12_5_5_5:5"
#その他の線種の例
# -W5t20_5:5
# -W5t15_5:5
# -W5t20_5_5_5_5_5:5
# -W5t30_5_5_5:5
# -W5t20_5_10_5:5
# 数字の出力フォーマット(project, grd2xzy等で使う)
gmtset D_FORMAT %lg #デフォルト(規定値)
# 桁数を指定（例：123.45678)
#gmtset D_FORMAT %12.6f
# 状況に応じて
#gmtset D_FORMAT %12.5g
----------------------
End of gmtpar.sh
----------------------
======================
note.sh
======================
# Print time, current working directory and output filename
export LANG=C
currentdir=`pwd`
if [ ${#currentdir} -gt 90 ]; then
curdir1=${currentdir:1:90}
curdir2=${currentdir:91}
else
curdir1=$currentdir
curdir2="\ "
fi
now=`date`
host=`hostname`
#comment=" "
time1=$(ls -l $in | awk '{print $6, $7, $8}')
pstext -JX6/1.2 -R0/1/0/1.2 -N -X${xoffset:-0} -Y${yoffset:-0} << EOF -O >> $out
0 1.1 9 0 1 LM $0 $@
0 0.95 9 0 1 LM ${now}
0 0.80 9 0 1 LM ${host}
0 0.65 9 0 1 LM ${curdir1}
0 0.50 9 0 1 LM Input: ${in} (${time1})
0 0.1 9 0 1 LM ${comment:-""}
EOF
# Output: ${out}
# 0 0.50 9 0 1 LM ${curdir2}
----------------------
End of note.sh
----------------------