c============================================ C This code reads sea level purresure dataset C from NOAA & cal Southern Oscillation Index! C K.Mizobata C 2023.10.15 c============================================ C234567-------------------------------------- character*120 header integer iyear(72) real Dar(12,72),Tah(12,72) real mslpd(12),mslpt(12) real aDar(12,72),aTah(12,72) real sDar,sTah real nDar(12,72),nTah(12,72) real MSD,SOI(12,72) real SOI2(12*72),date2(12*72) real SOI3(744),SOI121(744) ! for 121month moving ave real SOI4(732),SOI13(732) ! for 121month moving ave !-------------------------------------------- !-- read SLP(1000mb subtracted) @ Darwin ---- !-------------------------------------------- open(10,file="./Darwin.txt") read(10,'(a120)') header !<--ヘッダー空読み !write(*,'(a120)') header i=1 do j=1,72 ! read(10,*) iyear(j),Dar(i,j),Dar(i+1,j), ! & Dar(i+2,j),Dar(i+3,j),Dar(i+4,j),Dar(i+5,j), ! & Dar(i+6,j),Dar(i+7,j),Dar(i+8,j),Dar(i+9,j), ! & Dar(i+10,j),Dar(i+11,j) read(10,*) iyear(j),(Dar(i,j),i=1,12) 100 format(i4,12f6.1) enddo close(10) !-------------------------------------------- !-- read SLP(1000mb subtracted) @ Tahiti ---- !-------------------------------------------- open(11,file="./Tahiti.txt") read(11,'(a120)') header i=1 do j=1,72 read(11,*) iyear(j),(Tah(i,j),i=1,12) enddo close(11) !------------------------------------ ! calculate mean slp for each month ! CAUTION! PERIOD is strictly determined by NOAA ! FROM 1981 to 2010 !------------------------------------ do i=1,12 mslpt(i)=0; mslpd(i)=0 do j=1,72 if(iyear(j).ge.1981.and.iyear(j).le.2010)then mslpt(i)=mslpt(i)+Tah(i,j) ! sum mslpd(i)=mslpd(i)+Dar(i,j) endif enddo mslpt(i)=mslpt(i)/(2010-1981+1) !mean mslpd(i)=mslpd(i)/(2010-1981+1) ! write(*,*)'k=',k,'mslpd=',mslpd(k) ! write(*,*)'i=',i,'mslpt=',mslpt(i) enddo !------------------------------------ ! calculate SLP deviation from mean !------------------------------------ do i=1,12 do j=1,72 aTah(i,j)=Tah(i,j)-mslpt(i) aDar(i,j)=Dar(i,j)-mslpd(i) enddo enddo ! write(*,*)" " ! write(*,*) "SLP Anomaly @ Darwin" ! do j=1,72 ! write(*,100) iyear(j),(aDar(i,j),i=1,12) ! enddo ! write(*,*)" " ! write(*,*) "SLP Anomaly @ Tahiti" ! do j=1,72 ! write(*,100) iyear(j),(aTah(i,j),i=1,12) ! enddo !------------------------------------ ! calculate SLP STANDRAD DEV !------------------------------------ sTah=0; sDar=0; icount=0 do j=1,72 if(iyear(j).ge.1981.and.iyear(j).le.2010)then do i=1,12 sTah=sTah+aTah(i,j)**2 sDar=sDar+aDar(i,j)**2 icount=icount+1 enddo endif enddo sTah=sqrt(sTah/icount) sDar=sqrt(sDar/icount) write(*,*)" " write(*,*) "Standard dev. of SLP @ Tahiti=",sTah write(*,*) "Standard dev. of SLP @ Darwin=",sDar write(*,*) "icount=",icount !------------------------------------ ! Normalized SLP @ Darwin and Tahiti !------------------------------------ !wirte(*,*) "Normalized SLP@ Tahiti" do j=1,72 do i=1,12 !nTah(i,j)=(Tah(i,j)-mslpt(i))/sTah !nDar(i,j)=(Dar(i,j)-mslpd(i))/sDar nTah(i,j)=aTah(i,j)/sTah nDar(i,j)=aDar(i,j)/sDar enddo !write(*,100) iyear(j),(nTah(i,j),i=1,12) enddo !------------------------------------ ! Monthly mean standard deviation (MSD) !------------------------------------ MSD=0; icount2=0 do j=1,72 if(iyear(j).ge.1981.and.iyear(j).le.2010)then do i=1,12 MSD=MSD+(nTah(i,j)-nDar(i,j))**2 icount2=icount2+1 enddo endif enddo MSD=sqrt(MSD/real(icount2)) write(*,*)" " write(*,*)"MSD=",MSD write(*,*)"icount2=",icount2 !------------------------------------ ! Final stage! Cal SOI !------------------------------------ do j=1,72 do i=1,12 SOI(i,j)=(nTah(i,j)-nDar(i,j))/MSD enddo enddo !------------------------------------------- ! Timeseries SOI will be used later ! to investigate the relationship ! between SOI and SST@Nino monitoring areas !------------------------------------------- decim=1./24.!(365/12/2/365) decim2=decim*2 icount3=1 open(15,file="SOItimeseries.txt") do j=1,72 do i=1,12 k=i+(j-1)*12 SOI2(k)=SOI(i,j) date=real(iyear(j))+decim+decim2*(i-1) date2(k)=date !write(*,200) iyear(j),i,k,SOI2(k) write(15,200) date,SOI2(k) enddo enddo 200 format(2f12.4) close(15) !<-- caution: should be closed ! otherwise whole data cannot be saved ! !----------------------------------- ! Output !----------------------------------- write(*,*)" " write(*,*)"Southern Oscillation Index" open(20,file="./SOI.txt") do j=1,72 write(20,100) iyear(j),(SOI(i,j),i=1,12) ! write(*,110) iyear(j),(SOI(i,j),i=1,12) 110 format(i4,12f6.1) enddo close(20) !----------------------------------- ! 121-month moving average !----------------------------------- do k=ここは自分で考えましょう ! in reality,from 61 to 804 dummy=0 do i=1,121 dummy=dummy+SOI2(ここは自分で考えましょう) SOI121(k)=dummy/121 enddo enddo !--------------------------------- ! remove 121-month mov ave from original !--------------------------------- do k=ここは自分で考えましょう SOI3(k)=SOI2(ここは自分で考えましょう)-SOI121(k) enddo !----------------------------------- ! 13-month moving average !----------------------------------- do k=ここは自分で考えましょう ! in reality,from 67 to 798 dummy2=0 do i=1,13 dummy2=dummy2+SOI3(k+i-1) SOI13(k)=dummy2/13 enddo enddo !------ output smothing data ------! open(20,file="SOIsmoothed.txt") do k=ここは自分で考えましょう write(20,200) date2(ここは自分で考えましょう),SOI13(k) enddo close(20) !<-- caution: should be closed ! otherwise whole data cannot be saved ! stop !<---- we need to add it ALWAYS! end