コンポジット解析
（参考）
https://www.nature.com/articles/srep33790
上記リンク先は数値モデルのデータを使い、出力されたデータから"Strong El Nino”期および"Strong La Nina”期における海面温度を計算している。
ここでは同様の計算を観測値を用いて計算します。
（実践）
以前使った、エルニーニョ監視海域Nino3.4における水温データを指標にして、"Strong El Nino”期および"Strong La Nina”期における平均海面温度の計算（コンポジット解析）を行う
今回はCaiらがよく使う指標で「0.5 s.d.」を用いる。s.d.は標準偏差standard deviationのこと。
例えば
More extreme swings of the South Pacific convergence zone due to greenhouse warming | Nature
ここでは、Nino3.4海域の海面温度時系列の平均と標準偏差を用いた閾値を計算する。
用意するデータ（いずれも前回使用したもの）
・海面温度データOISST
・エルニーニョ監視海域のSSTアノマリー
c234567---------------------------------
character*120 ifl(494)
integer num
integer*1 lmask(360,180) ! land mask
real dummy(360,180),msst(360,180)
real sst3d(360,180,494),xvar(494),dnum(360,180)
real YR(502),MO(502)
real nino(502,8)
real ELNINO(360,180),LANINA(360,180)
real COUNT1(360,180),COUNT2(360,180)
! SST : Dec. 1981 to Jan. 2023 (2~494)
! Nino SST: Jan. 1982 to Oct. 2023 (1~493)
C---------------------------------------
C read the list of FILENAME
C and make time axis
C---------------------------------------
open(10,file="list.txt")
read(10,*) num
do k=1,num
read(10,'(a120)') ifl(k)
write(*,'(a120)') ifl(k)
xvar(k)=k
enddo
close(10)
C---------------------------------------
C read SST files and make 3-D dataset
C---------------------------------------
do k=1,num
open(15,file=ifl(k),form='unformatted',
& access="direct",status="old",recl=360*180*4)
read(15,rec=1) dummy
do j=1,180
do i=1,360
sst3d(i,j,k)=dummy(i,j)
enddo
enddo
close(15)
enddo
C---------------------
C read LAND MASK
C---------------------
open(20,file='lmask.bin',form='unformatted',
& access="direct",status="old",recl=360*180)
read(20,rec=1) lmask
close(20)
C----------------------------
C cal mean for anomaly
C----------------------------
do j=1,180
do i=1,360
msst(i,j)=0; ic=0
do k=1,num
if(sst3d(i,j,k).gt.-1.79)then
msst(i,j)=msst(i,j)+sst3d(i,j,k)
ic=ic+1
endif
enddo
msst(i,j)=msst(i,j)/ic
dnum(i,j)=real(ic)
enddo
enddo
!------- read SST timeseries ---------!
inum2=1
open(11,file="sstoi.indices.txt")
read(11,'(a120)') header
write(*,'(a120)') header
1001 read(11,*,end=666) YR(inum2),MO(inum2),(nino(inum2,k),k=1,8)
write(*,201) YR(inum2),MO(inum2),(nino(inum2,k),k=1,8)
inum2=inum2+1
goto 1001
666 inum2=inum2-1
close(11)
write(*,*) "number of data(SST) = ",inum2
201 format(f6.0,1x,f4.0,1x,8f8.2)
C==================================
C define "Strong" Nino using Nino3.4
C cal. standard deviation
C==================================
mNino34=0; ! mean vaue of Nino3.4 <--initialize!
do i=1,inum2
mNino34=mNino34+nino(i,8);
enddo
mNino34=mNino34/real(inum2);
sNino34=0; ! std of Nino3.4 <--initialize!
do i=1,inum2
sNino34=sNino34+(nino(i,8)-mNino34)**2;
enddo
sNino34=sNino34/real(inum2); ! <--- variance
sNino34=sqrt(sNino34); ! <--- std
cr1=mNino34+0.5*sNino34 !<--criteria for strong El Nino
cr2=mNino34-0.5*sNino34 !<--criteria for strong La Nina
write(*,*) cr1,cr2
C============================================
C COMPOSITE SST based on criteria above
C============================================
! SST : Dec. 1981 to Jan. 2023 (2~494)
! Nino SST: Jan. 1982 to Oct. 2023 (1~493)
do j=1,180
do i=1,360
ELNINO(i,j)=0; COUNT1(i,j)=0; ic=0
LANINA(i,j)=0; COUNT2(i,j)=0;
enddo
enddo
do k=2,num
if(nino(k-1,8).gt.cr1)then
write(*,*) "Strong k=",k
do j=1,180
do i=1,360
if(sst3d(i,j,k).gt.-1.79)then
ELNINO(i,j)=ELNINO(i,j)+sst3d(i,j,k)-msst(i,j)
COUNT1(i,j)=COUNT1(i,j)+1
endif
enddo
enddo
endif
if(nino(k-1,8).lt.cr2)then
do j=1,180
do i=1,360
if(sst3d(i,j,k).gt.-1.79)then
LANINA(i,j)=LANINA(i,j)+sst3d(i,j,k)-msst(i,j)
COUNT2(i,j)=COUNT2(i,j)+1
endif
enddo
enddo
endif
enddo
! finalize !
do j=1,180
do i=1,360
if(COUNT1(i,j).gt.0)then
ELNINO(i,j)=ELNINO(i,j)/COUNT1(i,j);
endif
if(COUNT2(i,j).gt.0)then
LANINA(i,j)=LANINA(i,j)/COUNT2(i,j);
endif
enddo
enddo
C----------------------------
C ERASE SST data on "LAND"
C----------------------------
do j=1,180
do i=1,360
if(lmask(i,j)==0)then
ELNINO(i,j)=0
LANINA(i,j)=0
endif
enddo
enddo
C----------------------------
C OUTPUT-->write BINARY file
C----------------------------
open(80,file='ELNINO.bin',form='unformatted',
& access="direct",recl=360*180*4)
write(80,rec=1) ELNINO
open(81,file='LANINA.bin',form='unformatted',
& access="direct",recl=360*180*4)
write(81,rec=1) LANINA
stop
end
gmt xyz2grd ELNINO.bin -R0.5/359.5/-89.5/89.5 -I1 -ZTLf -Gelnino -di-999
gmt makecpt -Chaxby -T-1/1/0.05 -N
gmt grdimage elnino -R0.5/359.5/-89.5/89.5 -JR16c
gmt grdcontour elnino -C0.5 -A1+f10p -Gd8c -S4
gmt coast -B+t"SST during strong El Nino" -B -Gwhite -Wthinnest
gmt colorbar -DjTR+o0c+w3c/0.4c -Bx -By+loC -I -F+gwhite+p1p
gmt end show