11:データ観察(時空間解析2) 各海域における海面高度の上昇傾向 今回はCMEMSが配布する海面高度データを使用する。 海面温度データと同じく、360x180の1度グリッドデータとしている。 以下の図は、季節変動を除去したのちに、線形トレンドを計算してプロットしている。 なお、1993年~2000年までの上昇トレンドは1度グリッドでは2.66mm/年となった。 2013年以降の上昇トレンドを計算してみよう。 似たような計算結果はここでもみられる。ただしグリッドサイズが異なるので結果は若干違う。 計算の手順 データはここ! 1)1993年1月~2023年5月までの各月における平均海面水位(実際には、1993年~2012年の平均場からのズレとなる海面偏差の平均)を計算する。→平均海面水位の時系列の取得 2)平均海面水位を従属変数、時間を説明変数として線形一次の回帰式を求める。 Y(海面水位)=a x X(時間, xvar)+b 3)時間xvarと、求めた傾きaと切片b、で予測される1993年~2023年までの海面水位を求める。 4)元の時系列に季節変動成分が含まれている。各月毎に、1)の時系列と3)の時系列の差の平均を求める。 サンプルコードでは変数monmeanとしている。 5)1)の時系列から4)を引くことで、季節変動シグナルを除去した海面水位時系列を得る。 6)5)の時系列に対して、任意の期間における線形一次の回帰式を求める。 得られた傾きはひと月での上昇率となるので、12を掛けて1年での上昇率に変換する。 |
サンプルコード c234567-------------------------------- character*23 ifl integer inum,inum2,smo,emo real dummy(360,180),msla real sla3d(360,180,365),xvar(365) real sla1d(360),sla1d2(365) real sla1df(365),sla1df2(365) real sla1d9303(120) real sla1d1323(125) real slope1,offset1 real slope2,offset2 real slope3,offset3 integer month(365),year(365) real monmean(12) !-------------------------------------- ! read SLA dataset; build SLA 3d(lon,lat,time) !-------------------------------------- inum=1 smo=1;emo=12; do iyr=1993,2023 if(iyr.eq.2023)then emo=5 endif do imo=smo,emo if(imo.ge.10)then write(ifl,'("./CMEMS/",i4,"_m",i2,"sla.bin")') iyr,imo else write(ifl,'("./CMEMS/",i4,"_m0",i1,"sla.bin")') iyr,imo endif write(*,200) ifl 200 format(a23) !--make file's name open(10,file=ifl,form='unformatted', & access="direct",status="old",recl=360*180*4) read(10,rec=1) dummy do j=1,180 do i=1,360 sla3d(i,j,inum)=dummy(i,j) enddo enddo year(inum)=iyr ! we will use it later! month(inum)=imo ! we will use it later! xvar(inum)=inum ! write(*,*) inum,year(inum),month(inum),xvar(inum) inum=inum+1 enddo enddo inum=inum-1 write(*,*) "inum=",inum !-------------------------------------- ! cal mean value of timeseries !-------------------------------------- do k=1,inum sla1d(k)=0; hoge=0 inum2=0; do j=1,180 !summation! do i=1,360 if(sla3d(i,j,k).gt.-900)then !sla1d(k)=sla1d(k)+sla3d(i,j,k) hoge=hoge+sla3d(i,j,k) inum2=inum2+1 endif enddo enddo sla1d(k)=hoge/real(inum2) ! average !write(*,*) k,inum2,sla1d(k) enddo !---------------------------- ! Least-square fitting (1-D) !---------------------------- B=0; C=0; D=0; E=0; do k=1,inum B=B+ここは考えましょう C=C+ここは考えましょう D=D+ここは考えましょう E=E+ここは考えましょう enddo slope=ここは考えましょう offset=ここは考えましょう C------------------------------ ! Build SLA timeseries based on results above ! sla1d2 is the linear trend of SLA C------------------------------ do k=1,inum sla1d2(k)=ここは考えましょう enddo !----------------------------------------------- ! Deseasonalized but keep the number of data! !----------------------------------------------- ! 1st step: cal seasonal variation of sla relative to linear trend ! do i=1,12 ! icount=0 ! monmean(i)=0 ! do k=1,inum ! if(ここは考えましょう)then ! monmean(i)=monmean(i)+(ここは考えましょう) ! icount=icount+1 ! endif ! enddo ! monmean(i)=monmean(i)/real(icount) ! enddo ! do i=1,12 ! write(*,*) "monmean(",i,")=",monmean(i) ! enddo ! 2nd step: erase seasonal trend relative to sla trend ! do i=1,12 ! do k=1,inum ! if(ここは考えましょう)then ! sla1df(k)=ここは考えましょう ! endif ! enddo ! enddo ! ! !---------------------------- ! Least-square fitting (1-D) applying to "sla1df" !---------------------------- ! B=0; C=0; D=0; E=0; ! do k=1,inum ! B=B+ここは考えましょう ! C=C+ここは考えましょう ! D=D+ここは考えましょう ! E=E+ここは考えましょう ! enddo ! slope1=ここは考えましょう ! offset1=ここは考えましょう ! write(*,*) "slope1,offset1",slope1*12,offset1 C------------------------------ ! Build SLA timeseries based on results above ! sla1df2 is the linear trend of SLA, WITHOUT SEASONAL VARIATION C------------------------------ ! do k=1,inum ! sla1df2(k)=ここは考えましょう ! enddo !---------------------------- C OUTPUT-->write ascii file C---------------------------- decim=1./24.!(365/12/2/365) decim2=decim*2 yr=1993 mo=1 open(80,file='GlobalmeanSLA.txt') open(81,file='GlobalmeanSLA2.txt') do k=1,inum date=yr+decim+decim2*(mo-1) write(80,100) date,sla1d(k)*100 write(81,100) date,sla1d2(k)*100 mo=mo+1 if(mo.eq.13)then mo=1 yr=yr+1 endif enddo close(80) close(81) 100 format(2f10.4) stop end |
地衡流、渦運動エネルギー、渦度の計算 参考資料 |