* PACKAGE DSM  !" SMPACK を利用した球面調和関数変換 ( 多層 )
*
*" [HIS] 98/01/18 (akahori) 新規作成 SMPACK-0.0 対応
*"       98/03/07 (akahori) SMPACK-1.0 対応
*
*=====================================================================
      SUBROUTINE DSMINI     !" 球関数(多層)を初期化 
      IMPLICIT NONE

*   [PARAM] 
#ifdef SYS_IBMS
      INCLUDE   (ZCDIM)                      !" 格子点数，波数
      INCLUDE   (ZCCOM)                      !" 標準物理定数
#else
#include        "zcdim.F"
#include        "zccom.F"
#endif

! [INTERNAL WORK FOR SMPACK]

      ReaL*8    X   (IMAX)                !" 南北格子点座標値(radian)
      ReaL*8    Y   (JMAX)                !" 東西格子点座標値(radian)
      ReaL*8    W   (JMAX)                !" ガウス荷重
      INTEGER   IT  (5)                   !" SMT???が用いる配列
      ReaL*8    T   (IMAX*2)              !" SMT???が用いる配列
      INTEGER   IP  (KMAX*(2*MMAX+1))     !" SMT???が用いる配列
      ReaL*8    P   (KMAX*(MMAX+1)*(MMAX+3)*JMAX/2)
                                          !" SMT???が用いる配列
      ReaL*8    QSG (JMAX/2)              !" SMTS2Gが用いる配列
      ReaL*8    QGS (JMAX/2)              !" SMTG2Sが用いる配列
      ReaL*8    QSV (JMAX/2)              !" SMTS?Vが用いる配列
      ReaL*8    QVS (JMAX/2)              !" SMTV?Sが用いる配列
      ReaL*8    QSU (JMAX/2)              !" SMTS?Vが用いる配列
      ReaL*8    QUS (JMAX/2)              !" SMTV?Sが用いる配列
      ReaL*8    R   (KMAX*MMAX*(MMAX-1)*2)!" SMTV?S,SMTS?Vが用いる配列
      INTEGER   ML  (NMDIM)               !" 東西波数 m のテーブル
                                          !" SMTV?S, SMTS?Vが用いる
      INTEGER   NL  (NMDIM)               !" 全波数 n のテーブル
      ReaL*8    RF  (NMDIM)               !" SMCLAPが用いる配列
      ReaL*8    RB  (NMDIM)               !" SMCLAPが用いる配列

      SAVE      X,Y,W,IT,T,IP,P
      SAVE      QSG,QGS,QSV,QVS,QSU,QUS
      SAVE      R,ML,NL,RF,RB

      ReaL*8    QSVMER(JMAX/2)            !" QSV*ER
      ReaL*8    QVSDER(JMAX/2)            !" QVS/ER
      ReaL*8    QSUMER(JMAX/2)            !" QSU*ER
      ReaL*8    QUSDER(JMAX/2)            !" QUS/ER
      SAVE      QSVMER, QVSDER, QSUMER, QUSDER

! [INPUT or OUTPUT]
      REAL      GDATA( IDIM,JDIM,KMAX )
      REAL      GDATB( IDIM,JDIM,KMAX )
      REAL      WDATA( NMDIM,KMAX )
      REAL      WDATB( NMDIM,KMAX )

! [INTERNAL WORK]
      INTEGER   J, K, NM
      CHARACTER HFUNC*4

      ReaL*8    WORKA( NMDIM,KMAX )       !" ワーク
      ReaL*8    WORKB( NMDIM,KMAX )       !" ワーク
      ReaL*8    WORK ( IDIM*JDIM*KMAX*3 ) !" ワーク
      ReaL*8    GW   ( JMAX )             !" GW ( ガウス荷重 )
      ReaL*8    Y0   ( JMAX )             !" 東西格子点座標値(radian)

      LOGICAL   OFIRST
      DATA      OFIRST /.TRUE./
      SAVE      OFIRST

!"  < このサブルーチンは一度だけ呼べばよい >

      IF ( .NOT. OFIRST ) RETURN
      WRITE ( 6,* ) ' SMPACK-1.0 FOR MULTI-LAYERS DATE=98/03/03'
      OFIRST = .FALSE.

!"  < SMPACK の使用条件をチェック >

      IF(MMAX.LT.1)      CALL MSGDMP('E','DSMINI', 'MMAX < 1' )
      IF(KMAX.LT.1)      CALL MSGDMP('E','DSMINI', 'KMAX < 1' )
      IF(JMAX.LE.MMAX)   CALL MSGDMP('E','DSMINI', 'JMAX <= MMAX' )
      IF(IMAX.LE.2*MMAX) CALL MSGDMP('E','DSMINI', 'MMAX <= 2*MMAX')

      IF(MINT.NE.1 )     CALL MSGDMP('E','DSMINI', 'MINT .NE. 1' )
      IF(MMAX.NE.LMAX )  CALL MSGDMP('E','DSMINI', 'MMAX .NE. LMAX')
      IF(JMAX.LT.2 )     CALL MSGDMP('E','DSMINI', 'JMAX < 2')
      IF(IMAX.LT.2 )     CALL MSGDMP('E','DSMINI', 'IMAX < 2')
      IF(NMDIM.NE.(MMAX+1)*(MMAX+1)) CALL MSGDMP('E','DSMINI','NMDIM')

!"  < 球面調和関数初期化 >

      CALL SMINIT
     I     (MMAX,IMAX,JMAX,KMAX,
     O      X,Y,W,IT,T,IP,P,
     O      QSG,QGS,QSV,QVS,QSU,QUS,
     O      R,ML,NL,RF,RB )

!"  < 地球半径ファクターを調整 >

      DO J=1,JMAX/2
        QSVMER(J)=QSV(J)*ER
        QVSDER(J)=QVS(J)/ER
        QSUMER(J)=QSU(J)*ER
        QUSDER(J)=QUS(J)/ER
      ENDDO

      RETURN

**********************************************************************
      ENTRY DSMG2S          !" 格子→スペクトル
     O         ( WDATA ,
     I           GDATA )

      CALL SMTG2S
     I     ( MMAX, IMAX, IDIM, JMAX, JDIM, KMAX,
     I       GDATA,
     O       WDATA, 
     W       WORK,
     I       IT, T, IP, P, QGS ) 

      RETURN

**********************************************************************
      ENTRY DSMV2S          !" 格子(ベクトル)→スペクトル(発散・回転)
     O         ( WDATA , WDATB,
     I           GDATA , GDATB,
     I           HFUNC  )

      IF ( HFUNC(1:1) .EQ. ' ' ) THEN
        CALL SMTV2S
     I     ( MMAX, IMAX, IDIM, JMAX, JDIM, KMAX,
     I       GDATA, GDATB,
     O       WDATA, WDATB,
     W       WORK,
     I       IT, T, IP, P, QVSDER, R, ML ) 
      ELSE
        CALL SMTV2S
     I     ( MMAX, IMAX, IDIM, JMAX, JDIM, KMAX,
     I       GDATA, GDATB,
     O       WDATA, WDATB,
     W       WORK,
     I       IT, T, IP, P, QUSDER, R, ML ) 
      ENDIF

      RETURN

**********************************************************************
      ENTRY DSMV1S          !" 格子(ベクトル)→スペクトル(発散)
     M         ( WDATA ,
     I           GDATA , GDATB,
     I           HFUNC  )

      IF ( HFUNC(4:4) .EQ. 'S' ) THEN
        CALL COPY ( WORKA, WDATA, NMDIM*KMAX )
      ENDIF

      IF ( HFUNC(1:1) .EQ. ' ' ) THEN   !" AGCM5 ではつかっていない
        CALL SMTV1S             
     I     ( MMAX, IMAX, IDIM, JMAX, JDIM, KMAX,
     I       GDATA, GDATB,
     O       WDATA,
     W       WORK,
     I       IT, T, IP, P, QVSDER, R, ML ) 
      ELSE
        CALL SMTV1S
     I     ( MMAX, IMAX, IDIM, JMAX, JDIM, KMAX,
     I       GDATA, GDATB,
     O       WDATA,
     W       WORK,
     I       IT, T, IP, P, QUSDER, R, ML ) 
      ENDIF

      IF ( HFUNC(4:4) .EQ. 'S' ) THEN
      DO K=1,KMAX
        DO NM=1,NMDIM
          WDATA(NM,K)=WORKA(NM,K) - WDATA(NM,K)
        ENDDO
      ENDDO
      ENDIF

      RETURN

**********************************************************************
      ENTRY DSMS2G          !" スペクトル→格子
     O         ( GDATA ,
     I           WDATA  )

      CALL SMTS2G
     I     ( MMAX, IMAX, IDIM, JMAX, JDIM, KMAX,
     I       WDATA,
     O       GDATA, 
     W       WORK,
     I       IT, T, IP, P, QSG ) 

      RETURN

**********************************************************************
      ENTRY DSMS2V          !" スペクトル(発散・回転)→格子(ベクトル)
     O         ( GDATA , GDATB,
     I           WDATA , WDATB,
     I           HFUNC         )

      CALL SMCLAP ( MMAX, KMAX, WDATA, WORKA, RB )
      CALL SMCLAP ( MMAX, KMAX, WDATB, WORKB, RB )

      IF ( HFUNC(1:1) .EQ. ' ' ) THEN

        CALL SMTS2V
     I       ( MMAX, IMAX, IDIM, JMAX, JDIM, KMAX,
     I         WORKA, WORKB,
     O         GDATA, GDATB,
     W         WORK,
     I         IT, T, IP, P, QSVMER, R, ML ) 

      ELSE                             !" AGCM5 ではつかっていない

        CALL SMTS2V
     I       ( MMAX, IMAX, IDIM, JMAX, JDIM, KMAX,
     I         WORKA, WORKB,
     O         GDATA, GDATB,
     W         WORK,
     I         IT, T, IP, P, QSUMER, R, ML ) 

      ENDIF

      RETURN

***************************************************************      
      ENTRY  DSMGAU          !" ガウス緯度・ガウス荷重の取得
     O         ( Y0 , GW ) 

      CALL COPY ( Y0, Y, JMAX ) 
      CALL COPY ( GW, W, JMAX ) 

      RETURN
      END









