2018-06-08_15-09
/work2/am/2017.LLC.Pacific/HAGIWARA_MASK/OISST.MONTH.AVE.STD
am@localhost
$ srcdump.sh OISST.MONTH.AVE.STD.sh OISST.MONTH.AVE.STD.f90 OISST.MONTH.AVE.STD.PL.sh >l.txt
#------------------------------
# List of the following files:
#------------------------------
OISST.MONTH.AVE.STD.sh
OISST.MONTH.AVE.STD.f90
OISST.MONTH.AVE.STD.PL.sh
#------------------------------
# Machine info
#------------------------------
localhost
/work2/am/2017.LLC.Pacific/HAGIWARA_MASK/OISST.MONTH.AVE.STD
Fri Jun 8 15:09:38 JST 2018
#======================
# OISST.MONTH.AVE.STD.sh
#======================
#!/bin/bash
src=$(basename $0 .sh).f90
if [ ! -f $src ]; then
echo
echo ERROR in $0 : NO SUCH FILE, $src
echo
exit 1
fi
#============================
yyyys=2016
mms=06
dds=05
yyyye=2016
mme=06
dde=25
wlon=140
elon=180
slat=30
nlat=50
#============================
outdate=${yyyy}${mm}${dd}
indir=/misc/raid112MANDA/KakenA.H28/OISST/
infle=${indir}/sst.day.mean.${yyyys}.v2.nc
if [ ! -f $infle ]; then
echo
echo ERROR in $0 : No such file, $infle
echo
echo exit
fi
nml=$(basename $0 .sh).nml
cat<<EOF>$nml
¶
infle="${infle}"
yyyys=${yyyys}
mms=${mms}
dds=${dds}
yyyye=${yyyye}
mme=${mme}
dde=${dde}
wlon=${wlon}
elon=${elon}
slat=${slat}
nlat=${nlat}
&end
EOF
echo
ls -lh $nml
echo
cat $nml
echo
exe=$(basename $0 .sh).exe
DEBUG="-CB -traceback -fpe0"
BINOUT="-convert big_endian -assume byterecl"
OPT="-L/usr/local/netcdf-4.1.3/lib -lnetcdf -lnetcdff -I/usr/local/netcdf-4.1.3/include"
ifort $src -o $exe $DEBUG $OPT
if [ $? -ne 0 ]; then
echo
echo COMPILE ERROR!
echo
exit 1
else
echo
echo COMPILE FINISHED.
ls -lh $exe
echo
fi
$exe <$nml
#----------------------
# End of OISST.MONTH.AVE.STD.sh
#----------------------
#======================
# OISST.MONTH.AVE.STD.f90
#======================
program main
! Description:
!
! Author: aym
!
! Host: aofd30
! Directory: /work1/aym/11.Work10/15.HeatBudjet/51.jcope2/32.NCEP_FLUX/11.test_read_flux/src
!
! Revision history:
! 2011-02-22 16:38
! Initial Version
! 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
! use
! implicit none
include '/usr/local/include/netcdf.inc'
! ncid:ファイルのID番号、 varid: 変 数のID番号
integer ncid,varid
character infle*500,ofle*500
character var*50
real,allocatable :: lon(:),lat(:)
real(8),allocatable :: time(:)
real,allocatable :: sdata(:,:,:)
! real,allocatable :: sdata(:,:,:)
real,allocatable :: buf(:,:,:)
real,allocatable :: sum(:,:),ave(:,:),sdv(:,:),trd(:,:),cnt(:,:)
real,allocatable :: P(:),Q(:) !FOR LINEAR REGRESSION
integer,allocatable :: dimids(:),nshape(:)
character*70 err_message
character*100 dimname !dimnameは各次元の名前。
real wlon,elon,slat,nlat
character(15)::domain !XXX-XXX_XXX_XXX
integer julday
integer jday0, jdays, jdaye
integer yyyys,mms,dds, yyyye,mme,dde, yyyyc,mmc,ddc
character(8)::sdate,edate
integer stat
namelist /para/infle,yyyys,mms,dds,yyyye,mme,dde,&
wlon,elon,slat,nlat
write(*,'(a)')'Program read_ncep starts.'
! write(*,*)''
read(*,nml=para)
jday0=julday(yyyys,01,01)
jdays=julday(yyyys,mms,dds)
jdaye=julday(yyyye,mme,dde)
! time index (n) for OISST dataset
ns=jdays-jday0+1
ne=jdaye-jday0+1
! CHECK
call caldat(ns+jday0-1,yyyyc,mmc,ddc)
print *,'ns=',ns
print *,yyyys,mms,dds, yyyyc,mmc,ddc
if(yyyys/=yyyyc .or. mms/=mmc .or. dds/=ddc)then
print *,'ERROR'
stop
endif
call caldat(ne+jday0-1,yyyyc,mmc,ddc)
print *,'ne=',ne
print *,yyyye,mme,dde, yyyyc,mmc,ddc
if(yyyye/=yyyyc .or. mme/=mmc .or. dde/=ddc)then
print *,'ERROR'
stop
endif
write(sdate(1:8),'(i4.4,i2.2,i2.2)')yyyys,mms,dds
write(edate(1:8),'(i4.4,i2.2,i2.2)')yyyye,mme,dde
im=1440
jm=720
nm=365
nn=ne-ns+1
allocate(lon(im),lat(jm),time(nm),sdata(im,jm,nm))
allocate(sum(im,jm),ave(im,jm),sdv(im,jm),trd(im,jm),buf(im,jm,nm))
allocate(cnt(im,jm))
print '(A,A)','Input: ',infle(1:lnblnk(infle))
stat = nf_open(infle, nf_nowrite, ncid) ! ファイルのopenとNetCDF ID(ncid)の取得
if(stat.ne.0) then
print *,' stat = nf_open(infle, nf_nowrite, ncid)'
print *,'ncid= ',ncid
print *,'status= ',stat
err_message = nf_strerror(stat)
write(6,*) err_message
stop
endif
var='sst'
! ファイルID(ncid)からvarで設定した変数のID(varid)を得る。
stat = nf_inq_varid(ncid,var,varid)
if(stat /= 0)then
print '(A)','stat = nf_inq_varid(ncid,var,varid)'
print '(A)',var(1:lnblnk(var))
print *,'varid=', varid
print *,'stat= ', stat
print *
stop
endif
! スケールファクターを得る。(*1)
! sta = nf_get_att_real(ncid,varid,'scale_factor',scale)
! print *,'scale= ',scale
! print *,'sta= ',sta
! print *
! オフセットの 値を得る(*2)
! sta = nf_get_att_real(ncid,varid,'add_offset',offset)
! print *,'offset=', offset
! print *,'sta= ',sta
! print *
! varidからndimsを得る。ndimsとは次元の数。二次元データなら2
stat = nf_inq_varndims(ncid, varid, ndims)
! varidからdimidsを得る。dimidsとはそれぞれの次元に対するID
allocate(dimids(ndims))
stat = nf_inq_vardimid(ncid, varid, dimids)
if(stat /= 0)then
print *,'stat = nf_inq_vardimid(ncid, varid, dimids)'
print *,'sta= ', sta
stop
endif
! do i=1,ndims
! write(6,'("dimids(",i2,")=",i9 )') i,dimids(i)
! enddo
! print *
allocate(nshape(ndims))
do i=1,ndims
! dimidsからnshapeを得る。
! nshapeとはそれぞれの次元に対する最大データ数 (格子数)
stat = nf_inq_dimlen(ncid, dimids(i), nshape(i))
! if(stat /= 0)then
write(6,'("nshape(",i2,")=",i9 )') i,nshape(i)
! stop
! endif
enddo !i
print *
do i=1,ndims
stat = nf_inq_dimname(ncid, dimids(i), dimname)
! if(stat /= 0)then
write(6,'("dimname(",i2,")=",1x,a23)') i,dimname
! endif
enddo
print *
print '(A,A)','Read ',var(1:lnblnk(var))
stat = nf_get_var_real(ncid, varid, sdata)
if(stat /= 0)then
print *,'stat= ',stat
stop
endif
var='lon'
print '(A,A)','READ ',var(1:lnblnk(var))
! ファイルID(ncid)からvarで設定した変数のID(varid)を得る。
stat = nf_inq_varid(ncid,var,varid)
if(stat /= 0)then
print *,'varid=', varid
print *,'sta= ', sta
print *,'Read ',var(1:lnblnk(var))
endif
stat = nf_get_var_real(ncid, varid, lon)
if(stat /= 0)then
print *,'stat= ',stat
endif
var='lat'
print '(A,A)','READ ',var(1:lnblnk(var))
stat = nf_inq_varid(ncid,var,varid)
if(stat /= 0)then
print *,'varid=', varid
print *,'sta= ', sta
print *,'Read ',var(1:lnblnk(var))
endif
stat = nf_get_var_real(ncid, varid, lat)
if(stat /= 0)then
print *,'stat= ',stat
print *
endif
buf(:,:,:)=0.0
sum(:,:) =0.0
ave(:,:) =0.0
sdv(:,:) =0.0
trd(:,:) =0.0
cnt(:,:) =0.0
print '(A)','SUM'
TIME_LOOP: do n=ns,ne
LAT_LOOP: do j=1,jm
LON_LOOP: do i=1,im
if(lon(i)>=wlon.and.lon(i)<=elon.and.lat(j)>=slat.and.lat(j)<=nlat)then
if(sdata(i,j,n)>0.0)then
buf(i,j,n)=sdata(i,j,n)
sum(i,j) =sum(i,j)+sdata(i,j,n)
cnt(i,j) =cnt(i,j)+1.0
endif !sdata
endif !lon lat
end do LON_LOOP
end do LAT_LOOP
end do TIME_LOOP
print '(A)','AVE'
do j=1,jm
do i=1,im
if(cnt(i,j)>10.0)then
ave(i,j)=sum(i,j)/cnt(i,j)
else
ave(i,j)=-999.
endif
end do !i
end do !j
print '(A)','SDV'
do j=1,jm
do i=1,im
if(cnt(i,j)>10.0)then
do n=ns,ne
if(buf(i,j,n)>0.0)then
sdv(i,j)=sdv(i,j)+(buf(i,j,n)-ave(i,j))**2
end if
end do !n
sdv(i,j)=sqrt(sdv(i,j)/(cnt(i,j)-1.0))
else
sdv(i,j)=-999.
end if
end do !i
end do !j
print '(A)','TREND'
do j=1,jm
do i=1,im
if(cnt(i,j)>10.0)then !cnt
ND=0
do n=ns,ne
if(buf(i,j,n)>0.0)then !buf
ND=ND+1
endif !buf
enddo !n
allocate(P(ND),Q(ND))
ntemp=0
do n=ns,ne !n
if(buf(i,j,n)>0.0)then !buf
ntemp=ntemp+1
P(ntemp)=float(ntemp)
Q(ntemp)=buf(i,j,n)
endif !buf
enddo !n
CALL LSQM(P,Q,ND,A,B)
trd(i,j)=B
deallocate(P,Q)
else !cnt
trd(i,j)=-999.
end if !cnt
end do !i
end do !j
write(domain( 1: 4),'(i3.3,A)')int(wlon),'-'
write(domain( 5: 8),'(i3.3,A)')int(elon),'_'
write(domain( 9:12),'(i3.3,A)')int(slat),'-'
write(domain(13:15),'(i3.3)')int(nlat)
ofle='OISST.MONTH.AVE.STD_'//trim(sdate)//'-'//trim(edate)//'_'//trim(domain)//'.txt'
open(20,file=ofle)
write(20,'(A,A)'),'# Input: ',infle(1:lnblnk(infle))
write(20,'(A,i4.4,1x,i2.2,1x,i2.2)'),'# yyyys mms dds= ',yyyys,mms,dds
write(20,'(A,i4.4,1x,i2.2,1x,i2.2)'),'# yyyye mme dde= ',yyyye,mme,dde
write(20,'(A,f10.3)'),'# wlon= ',wlon
write(20,'(A,f10.3)'),'# elon= ',elon
write(20,'(A,f10.3)'),'# slat= ',slat
write(20,'(A,f10.3)'),'# nlat= ',nlat
write(20,'(A,A)')'# lon(i),lat(j),ave(i,j),sdv(i,j),trd(i,j),cnt(i,j)'
do j=1,jm
do i=1,im
if(lon(i)>=wlon.and.lon(i)<=elon.and.lat(j)>=slat.and.lat(j)<=nlat)then
write(20,'(f9.3,f8.3,1x,3f10.3,f9.0,f10.1)') &
lon(i),lat(j),ave(i,j),sdv(i,j),trd(i,j),cnt(i,j)
endif
enddo !i
enddo !j
close(20)
print *
print '(A,A)','OUTPUT: ',trim(ofle)
print *
write(*,'(a)')'Done program'
! print '(A,A)','Output: ',ofle(1:lnblnk(ofle))
end program
! *1： scaleとは、ある変数の平均的なオーダー
! *2： offsetとは、ある変数の平均的な値
! data(i,j) = sdata(i,j) * scale + offset
! のようにして、本当のデータを得る
!caldat.f90
!ユリウス日→西暦
!
subroutine caldat(JULIAN,IYYY,MONTH,DD)
! Description:
!
! Author: am
!
! Host: aofd30
! Directory: /work2/am/12.Work11/75.Mitsui/52.Vertical_Profile/src
!
! Revision history:
! This file is created by /usr/local/mybin/nff.sh at 11:00 on 11-01-2011.
! use
! implicit none
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
!julday.f90
!西暦→ユリウス日
!
INTEGER FUNCTION JULDAY(IYYY,MONTH,DD)
INTEGER,intent(inout):: IYYY,MONTH,DD
!
! NAME IN/OUT DESCRIPTION
!
! IYYY I YEAR
! MONTH I MONTH (1 TO 12)
! DD I DAY OF MONTH
!
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
!***********************************************************************
!* LSQM.F
!* F:\0.Computer\TOOLBOX\BASIC_STATISTICS\回帰直線
!* CODED BY A.MANDA 1997.8
!* TASK
!C Linear Regression
!* PARAMETERS
!* (1)P : INPUT DATA [X-COOD.] (R)(INPUT)
!* (2)Q : INPUT DATA [Y-COOD.] (R)(INPUT)
!* (3)N : TOTAL DATA NUMBER (I)(INPUT)
!* (4)A : Y-SECTION OF LINE (R)(OUTPUT)
!* (5)B : GRADIENT OF LINE (R)(OUTPUT)
!* METHOD
!* SLAVE SUBROUTINE(S)
!* REMARK
!* REFERENCE(S)
!***********************************************************************
SUBROUTINE LSQM(P,Q,N,A,B)
!* INPUT
integer,intent(in)::N
REAL,intent(in):: P(N),Q(N)
!* OUTPUT
REAL,intent(inout):: A,B
TP=0.0
TQ=0.0
DO 10 I=1,N
TP=TP+P(I)
TQ=TQ+Q(I)
10 CONTINUE
PM=TP/FLOAT(N)
QM=TQ/FLOAT(N)
TPQ=0.0
TPP=0.0
DO 20 I=1,N
TPQ=TPQ+(P(I)-PM)*(Q(I)-QM)
TPP=TPP+(P(I)-PM)**2
20 CONTINUE
B=TPQ/TPP
A=QM-B*PM
RETURN
END SUBROUTINE LSQM
#----------------------
# End of OISST.MONTH.AVE.STD.f90
#----------------------
#======================
# OISST.MONTH.AVE.STD.PL.sh
#======================
#!/bin/bash
# Description:
#
# Author: am
#
# Host: localhost
# Directory: /work2/am/2017.KYUSHU-HOKUBU.HEAVY.RAIN/OISST.ASCII.OUT
#
# Revision history:
# This file is created by /usr/local/mybin/ngmt.sh at 09:13 on 02-12-2018.
. ./gmtpar.sh
echo "Bash script $0 starts."
gmtset HEADER_OFFSET = -0.3 HEADER_FONT_SIZE = 16p
range=140/180/30/50
size=Q130/5
xanot=a10f5
yanot=a5f5
in=OISST.MONTH.AVE.STD_20160605-20160625_140-180_030-050.txt
if [ ! -f $in ]; then
echo Error in $0 : No such file, $in
exit 1
fi
figdir= #"../fig/"
#if [ ! -d ${figdir} ];then
# mkdir -p $figdir
#fi
out=${figdir}$(basename $in .txt).ps
xoffset=1
yoffset=9.5
export LANG=C
curdir1=$(pwd)
now=$(date)
host=$(hostname)
time=$(ls -l ${in} | awk '{print $6, $7, $8}')
#time1=$(ls -l ${in1} | awk '{print $6, $7, $8}')
pstext <<EOF -JX6/1.5 -R0/1/0/1.5 -N -X${xoffset:-0} -Y${yoffset:-0} -K -P > $out
0 1.50 9 0 1 LM $0 $@
0 1.35 9 0 1 LM ${now}
0 1.20 9 0 1 LM ${host}
0 1.05 9 0 1 LM ${curdir1}
0 0.90 9 0 1 LM Input: ${in} (${time})
EOF
#0 0.75 9 0 1 LM Input:
#0 0.60 9 0 1 LM Input:
#0 0.45 9 0 1 LM Input:
title=AVE${sp}[${deg}C]
gmtset HEADER_OFFSET = 0. HEADER_FONT_SIZE = 16p
grd=$(basename $in .txt).AVE.grd
awk '{if($1!="#"&&$3>0)print $1,$2,$3}' $in|\
surface -R$range -I0.25/0.25 -T1 -G$grd
grdcontour $grd -R$range -J$size -W1 -C1 -Y-2.5 -O -K>>$out
grdcontour $grd -R$range -J$size -W3 -A5 -O -K>>$out
pscoast -R$range -J$size -W1 -G200 -O -K >>$out
anot=${xanot}/${yanot}:."${title}":WsNe
psbasemap -R -J$size -B$anot -O -K>>$out
title=SD${sp}[K]
gmtset HEADER_OFFSET = -0.3 HEADER_FONT_SIZE = 16p
grd2=$(basename $in .txt).SD.grd
cpt2=$(basename $in .txt).SD.cpt
makecpt -Chot -I -T0.5/3/0.5 >$cpt2
awk '{if($1!="#"&&$4>-999.)print $1,$2,$4}' $in|\
surface -R$range -I0.1/0.1 -T0.8 -G$grd2
grdimage $grd2 -R$range -J$size -C$cpt2 -Y-3 -O -K >>$out
grdcontour $grd -R$range -J$size -W6/${white} -C1 -O -K>>$out
grdcontour $grd -R$range -J$size -W2 -C1 -O -K>>$out
grdcontour $grd -R$range -J$size -W4 -A5 -O -K>>$out
psscale -D5.2/1.25/2.5/0.1 -C$cpt2 -E -Ba0.5f0.5 -O -K >>$out
pscoast -R$range -J$size -W1 -G200 -O -K >>$out
anot=${xanot}/${yanot}:."${title}":Wsne
psbasemap -R -J$size -B$anot -O -K>>$out
title=TREND${sp}[K/d]
grd2=$(basename $in .txt).TREND.grd
cpt2=$(basename $in .txt).TREND.cpt
makecpt -Cpolar -T-0.5/0.5/0.1 >$cpt2
awk '{if($1!="#"&&$5>-999.)print $1,$2,$5}' $in|\
surface -R$range -I0.1/0.1 -T0.8 -G$grd2
grdimage $grd2 -R$range -J$size -C$cpt2 -Y-3 -O -K >>$out
grdcontour $grd -R$range -J$size -W6/${white} -C1 -O -K>>$out
grdcontour $grd -R$range -J$size -W2 -C1 -O -K>>$out
grdcontour $grd -R$range -J$size -W4 -A5 -O -K>>$out
psscale -D5.2/1.25/2.5/0.1 -C$cpt2 -E -Ba0.5f0.1 -O -K >>$out
pscoast -R$range -J$size -W1 -G200 -O -K >>$out
anot=${xanot}/${yanot}:."${title}":Wsne
psbasemap -R -J$size -B$anot -O -K>>$out
rm -fv $(basename $in .txt)*.grd
#$(basename $in .txt)*.cpt
echo
echo Input: $in
echo
echo Output: $out
echo
echo "Done $0"
#----------------------
# End of OISST.MONTH.AVE.STD.PL.sh
#----------------------