!-----------------------------------------------------------------------
! 回転球面上のバロトロピック渦度方程式の時間発展     
!                                                2020/08/02 by K.Ishioka
!-----------------------------------------------------------------------
PROGRAM MODEL
  USE ISO_C_BINDING  
  IMPLICIT NONE
  INTEGER(8),PARAMETER :: IM=256, JM=128 ! λ, φ 方向の格子点数
  INTEGER(8),PARAMETER :: MM=85 ! 水平方向切断波数(微分とか除く)
  INTEGER(8),PARAMETER :: LM=(MM+1)*(MM+1) ! スペクトルを格納する配列の大きさ
  REAL(8) :: OMEGA, DT, DNU, TIM
  INTEGER(8) :: I, ITM, LEV, NDV
  REAL(8) :: VAR(LM),WORK(LM*3)
  REAL(8),ALLOCATABLE :: S(:),S1(:),S2(:)
  REAL(8) :: DRN(LM)
  REAL(8),PARAMETER :: PI=3.1415926535897932385D0 ! πの値
  !--- ISPACKの中で使われる配列
  INTEGER(8) :: IT(IM/2)  
  REAL(8) :: T(IM*3/2)
  REAL(8) :: R(((MM+1)*(2*(MM+1)-MM-1)+1)/4*3+(2*(MM+1)-MM)*(MM+1)/2+MM+1)
  INTEGER(8) :: JC(MM*(2*(MM+1)-MM-1)/16+MM)
  REAL(8) :: C(MM+1+MM*(2*MM-MM+1))
  REAL(8) :: D((MM+1+MM*(2*MM-MM+1)),2)
  REAL(8) :: YG(JM) ! ガウス緯度を格納する配列
  REAL(8),DIMENSION(:,:),POINTER:: P
  REAL(8),DIMENSION(:,:),POINTER:: G,G1,G2
  REAL(8),DIMENSION(:),POINTER:: W    
  TYPE(C_PTR) :: PP  
  TYPE(C_PTR) :: PG,PG1,PG2
  TYPE(C_PTR) :: PW

  !--- 配列の割り付け
  ALLOCATE(S(MM*(MM+4)+2))  ! S, S1, S2は緯度微分演算などのために大きめに
  ALLOCATE(S1(MM*(MM+4)+2)) ! とっていることに注意.
  ALLOCATE(S2(MM*(MM+4)+2)) !
  CALL MXALLC(PP,JM/2*(2*MM+5))
  CALL MXALLC(PG,JM*IM)
  CALL MXALLC(PG1,JM*IM)
  CALL MXALLC(PG2,JM*IM)    
  CALL MXALLC(PW,JM*IM*2)
  CALL C_F_POINTER(PP, P, [JM/2,2*MM+5])  
  CALL C_F_POINTER(PG, G, [IM,JM])
  CALL C_F_POINTER(PG1, G1, [IM,JM])
  CALL C_F_POINTER(PG2, G2, [IM,JM])    
  CALL C_F_POINTER(PW, W, [IM*JM*2])

  !---- 実験パラメターの設定 ---------------------------------------------
  
  OMEGA=2*PI  !自転角速度Ω
  ITM=100  !時間発展するステップ数
  NDV=20   ! Runge-Kuttaでのステップ分割数
  DT=0.2D0   ! ファイル出力の時間間隔
  LEV=3   !高階粘性項のラプラシアンの階数  
  DNU=10D0/(1D0*(MM*(MM+1)-2))**LEV !高階粘性係数の設定

  !---- モデルの初期化と従属変数の初期条件の設定 -------------------------
  CALL MODEL_INIT
  CALL VAR_INIT(VAR)

  !---- Runge-Kutta法による時間発展とファイルへの書き出し ----------------
  I=0
  TIM=0
  CALL ENEENS(I,VAR)
  OPEN(1,FILE='data.dat',FORM='UNFORMATTED')
  WRITE(1) VAR
  DO I=1,ITM
     CALL TXRKNU(LM,NDV,DT,TIM,VAR,WORK,TIME_DEVL,TIME_DERIVN)
     CALL ENEENS(I,VAR)
     WRITE(1) VAR     
  END DO
  CLOSE(1)

  !---  割り付けられた配列の解放
  DEALLOCATE(S)
  DEALLOCATE(S1)
  DEALLOCATE(S2)
  CALL MXFREE(PG)
  CALL MXFREE(PG1)
  CALL MXFREE(PG2)
  CALL MXFREE(PW)

CONTAINS    
!----------------------------------------------------------------------
! 時間発展に必要な配列等の初期化
!----------------------------------------------------------------------
  SUBROUTINE MODEL_INIT

    IMPLICIT NONE    
    INTEGER(8) :: L, J
    REAL(8) :: DTDNU

    CALL SXINI1(MM,MM+1,IM,IT,T,R)
    CALL SXINI2(MM,MM+1,JM,1_8,P,R,JC)    
    CALL SXINIC(MM,MM,C)
    CALL SXINID(MM,MM,D)

    ! ガウス緯度の情報の格納(φ)
    DO J=1,JM/2 
       YG(JM/2+J)=ASIN(P(J,1))
       YG(JM/2-J+1)=-ASIN(P(J,1))
    END DO

    DTDNU=0.5D0*DT/NDV*DNU
    DRN(1)=1
    DO L=1,LM
       DRN(L)=EXP(-DTDNU*(ABS(D(L,1))-2)**LEV) ! 高階粘性のための設定
    END DO

  END SUBROUTINE MODEL_INIT
!-----------------------------------------------------------------------
! 非線形項による時間変化率の計算
!-----------------------------------------------------------------------
  SUBROUTINE TIME_DERIVN(TIME,VAR,DVAR)

    IMPLICIT NONE        
    INTEGER(8) :: L
    REAL(8) :: TIME
    REAL(8) :: VAR(LM), DVAR(LM)

    CALL SXCLAP(MM,MM,VAR,S,D,2_8)
    CALL SXCS2X(MM,MM,S,S2)
    CALL SXCRPK(MM,MM,MM+1,S2,S1)    
    CALL SXCS2Y(MM,MM,S,S2,C)
    CALL SXTS2V(MM,MM+1,MM+1,IM,JM,S1,S2,G1,G2,IT,T,P,R,JC,W,1_8)
    ! G1: v, G2: -u
    
    DO L=1,LM
       S(L)=VAR(L)
    END DO
    S(2)=S(2)+2*OMEGA/SQRT(3D0) ! 渦度に惑星渦度の分を加える
    CALL SXTS2G(MM,MM+1,MM,IM,JM,S,G,IT,T,P,R,JC,W,0_8)
    ! G: q
    
    G1=G*G1 ! vq
    G2=G*G2 ! -uq
    CALL SXTV2S(MM,MM+1,MM+1,IM,JM,S1,S2,G1,G2,IT,T,P,R,JC,W,1_8)
    CALL SXCY2S(MM,MM,S1,S,C) ! S: ∂(cosφ vq)/cosφ∂φ
    CALL SXCRPK(MM,MM+1,MM,S2,S1)
    CALL SXCS2X(MM,MM,S1,S2) ! S2: ∂(-uq)/cosφ∂λ

    DO L=1,LM
       DVAR(L)=S2(L)-S(L)
    END DO

  END SUBROUTINE TIME_DERIVN
!-----------------------------------------------------------------------
! 線形項のみによる時間変化の計算
! (TIME から TIME + DTIME での VARの変化を与える. ここでは, TIMEと 
!  DTIME は実際には使われない)  
!-----------------------------------------------------------------------
  SUBROUTINE TIME_DEVL(TIME,DTIME,VAR)

    IMPLICIT NONE        
    INTEGER(8) :: L
    REAL(8) :: VAR(LM),TIME,DTIME

    DO L=1,LM
       VAR(L)=VAR(L)*DRN(L) ! 高階粘性による減衰
    END DO

  END SUBROUTINE TIME_DEVL
!-----------------------------------------------------------------------  
!     渦度場のスペクトルの初期値の設定
!-----------------------------------------------------------------------
  SUBROUTINE VAR_INIT(VAR)

    IMPLICIT NONE        
    INTEGER(8) :: L, I, J
    REAL(8) :: DA,DB,X0,Y0,CA,CB,CC,SA,SB,SC,CR,XG
    REAL(8) :: VAR(LM)

    DA=1D-3     !擾乱の振幅(α)
    DB=100D0    !擾乱の幅の逆数のルートに相当するパラメター(β)
    X0=0        !擾乱の中心位置のλ座標
    Y0=PI/4     !擾乱の中心位置のφ座標

    VAR=0
    VAR(4)=12*0.07D0*OMEGA ! 渦度の P^0_3 成分に 12*0.07D0 を与える

    ! 経度 X0, 緯度 Y0 を中心とするガウシアン的擾乱の設定
    CB=SIN(Y0)
    SB=COS(Y0)
    DO I=1,IM
       XG=2*PI*(I-1)/IM
       CA=COS(XG-X0)
       DO J=1,JM
          CC=SIN(YG(J))
          SC=COS(YG(J))
          CR=CB*CC+SB*SC*CA
          G(I,J)=DA*(EXP(DB*(CR-1))-(1-EXP(-2*DB))/(2*DB))
       END DO
    END DO
    CALL SXTG2S(MM,MM+1,MM,IM,JM,S,G,IT,T,P,R,JC,W,0_8)

    DO L=1,LM
       VAR(L)=VAR(L)+S(L) ! 擾乱を加える
    END DO
    
  END SUBROUTINE VAR_INIT
  !-----------------------------------------------------------------------
  !     全エネルギー, 全エンストロフィーの値の出力
  !-----------------------------------------------------------------------
  SUBROUTINE ENEENS(I,VAR)

    REAL(8) :: ENE, ENS
    REAL(8) :: VAR(LM)
    INTEGER(8) :: I, L

    ENE=0
    ENS=0
    DO L=2,MM+1
       ENE=ENE-D(L,2)*VAR(L)*VAR(L)
       ENS=ENS+VAR(L)*VAR(L)
    END DO
    DO L=MM+2,LM
       ENE=ENE-2*D(L,2)*VAR(L)*VAR(L)
       ENS=ENS+2*VAR(L)*VAR(L)
    END DO
    ENE=ENE/2
    ENS=ENS/2

    WRITE(6,'(A,I3,A,F21.15,A,F21.15)') &
         'step=',I,' energy=',ENE,' enstrophy=',ENS

  END SUBROUTINE ENEENS
  !--------------------------------------    
END program MODEL
!---------------------------------------------------------------------
! 非線形項と線形項を分離して Runge-Kutta法を使う
! (ISPACK1の TDRKNU の整数型変数のところを INTEGER(8)に変えたもの)
!---------------------------------------------------------------------
SUBROUTINE TXRKNU(N,M,H,X,Y,W,SUBL,SUBN)

  IMPLICIT NONE
  INTEGER(8) :: N,M,I
  REAL(8) :: H,X,Y(N),W(N,3),DX
  EXTERNAL SUBL,SUBN

  DX=H/M
  DO I=1,M
     CALL TXRKNL(N,DX,X,Y,W,SUBL,SUBN)
  END DO

END SUBROUTINE TXRKNU
!--------------------------------------        
SUBROUTINE TXRKNL(N,DX,X,Y,W,SUBL,SUBN)

  IMPLICIT NONE
  INTEGER(8) :: N, I
  REAL(8) :: DX,X,Y(N),W(N,3)
  EXTERNAL SUBL,SUBN

  CALL SUBN(X,Y,W(1,1))

  DO I=1,N
     W(I,2)=Y(I)+(DX/2)*W(I,1)
     W(I,3)=Y(I)+(DX/6)*W(I,1)
  END DO

  CALL SUBL(X,DX/2,Y)
  CALL SUBL(X,DX/2,W(1,2))
  CALL SUBL(X,DX/2,W(1,3))
  CALL SUBN(X+DX/2,W(1,2),W(1,1))

  DO I=1,N
     W(I,2)=Y(I)+(DX/2)*W(I,1)
     W(I,3)=W(I,3)+(DX/3)*W(I,1)
  END DO

  CALL SUBN(X+DX/2,W(1,2),W(1,1))

  DO I=1,N
     W(I,2)=Y(I)+DX*W(I,1)
     W(I,3)=W(I,3)+(DX/3)*W(I,1)
  END DO

  CALL SUBL(X+DX/2,DX/2,W(1,2))
  CALL SUBL(X+DX/2,DX/2,W(1,3))
  CALL SUBN(X+DX,W(1,2),W(1,1))

  X=X+DX
  DO I=1,N
     Y(I)=W(I,3)+(DX/6)*W(I,1)
  END DO

END SUBROUTINE TXRKNL
