* PACKAGE DDIAG  
*
*"  [HIS]  90/05/19(numaguti)
*
***********************************************************************
      SUBROUTINE DDIAG      
     I         ( GDU   , GDV   , GDT   , GDPS  , GDQ   ,
     I           GDVOR , GDDIV                          )
*





* PACKAGE  ZXT10L32 !" 格子点数，波数 X-Z
*
      INTEGER     IMAX           !" 経度格子数
      INTEGER     JMAX           !" 緯度格子数
      INTEGER     IDIM           !" 経度配列サイズ
      INTEGER     JDIM           !" 緯度配列サイズ
      INTEGER     MMAX           !" 最大東西波数
      INTEGER     LMAX           !" 最大南北波数
      INTEGER     NMAX           !" 最大全波数
      INTEGER     MINT           !" 東西波数間隔
      INTEGER     KMAX           !" 鉛直層数
      INTEGER     KDIM           !" 鉛直配列サイズ
      INTEGER     NMDIM          !" 水平波数成分数
      INTEGER     JMXHF          !" JMAX/2+1
      INTEGER     NWORK          !" ワーク領域大きさ
*
      PARAMETER ( IMAX  = 32,     JMAX =  1,
     &            IDIM  = IMAX+1, JDIM = JMAX,
     &            MMAX  = 10,     LMAX =  0,    NMAX = MMAX, MINT = 1,
     &            KMAX  = 32,     KDIM = KMAX,
     &            NMDIM = (MMAX/MINT+1)*(2*(NMAX+1)-MMAX)
     &                  - (NMAX-LMAX)/MINT*(NMAX-LMAX+1),
     &            JMXHF = JMAX/2+1,
     &            NWORK = IDIM*JDIM*KMAX*12                       )

* PACKAGE ZHDIM  !" 標準物理定数(common)  
*
*"  [HIS]  90/08/31(numaguti)
*
      ReaL*8       ER             !" 地球半径
      ReaL*8       GRAV                        !" 重力加速度
      ReaL*8       CP                          !" 大気定圧比熱
      ReaL*8       RAIR                        !" 大気気体定数
      ReaL*8       EL                          !" 水の凝結の潜熱
      ReaL*8       CPVAP                       !" 水蒸気定圧比熱
      ReaL*8       RVAP                        !" 水蒸気気体定数
      ReaL*8       DH2O                        !" 水の密度
      ReaL*8       EPSV                        !" 水と大気の分子量比ε
      ReaL*8       ES0                         !" ０℃飽和蒸気圧：Ｐａ
      ReaL*8       STB                         !" Stefan-Boltzman
      ReaL*8       FKARM                       !" Karman 定数
      ReaL*8       EPSVT                       !" 1/ε-1
*
      COMMON   /COMCON/
     &          ER    , GRAV  ,
     &          CP    , RAIR  ,
     &          EL    , CPVAP , RVAP  , DH2O  , EPSV  , ES0   ,
     &          STB   , FKARM , EPSVT

*
*   [INPUT]
      ReaL*8       GDU   ( IDIM, JDIM, KMAX )  
      ReaL*8       GDV   ( IDIM, JDIM, KMAX )  
      ReaL*8       GDT   ( IDIM, JDIM, KMAX )  
      ReaL*8       GDPS  ( IDIM, JDIM )        
      ReaL*8       GDQ   ( IDIM, JDIM, KMAX )  
      ReaL*8       GDVOR ( IDIM, JDIM, KMAX )  
      ReaL*8       GDDIV ( IDIM, JDIM, KMAX )  
*
*   [INTERNAL WORK]
      COMMON    /COMWRK /
     &           GDKE  , GDIE  , GDPE , GDLE ,
     &           GDTE  , GDENS , GDTM
*
      ReaL*8       GDKE  ( IDIM, JDIM, KMAX )  
      ReaL*8       GDIE  ( IDIM, JDIM, KMAX )  
      ReaL*8       GDPE  ( IDIM, JDIM, KMAX )  
      ReaL*8       GDLE  ( IDIM, JDIM, KMAX )  
      ReaL*8       GDTE  ( IDIM, JDIM, KMAX )  
      ReaL*8       GDENS ( IDIM, JDIM, KMAX )  
      ReaL*8       GDTM  ( IDIM, JDIM )        
*
      INTEGER    I,J,K

*    [INTERNAL ONCE]      
      LOGICAL    OFIRST
      DATA       OFIRST / .TRUE. /
      SAVE       OFIRST
*[ONCE]
      IF ( OFIRST ) THEN
         OFIRST = .FALSE.
         CALL DDIHRG
      ENDIF
*
*"         < 1. エネルギー，エンストロフィー >
*
*
      CALL GHYDRO
     O   ( GDPE  ,
     I     GDT     )
*
      DO 1100 K = 1, KMAX
         DO 1100 I = 1, IMAX
            DO 1100 J = 1, JMAX
*
               GDKE ( I,J,K )
     &                = ( GDU( I,J,K )**2 + GDV( I,J,K )**2 )/ 2.  
     &                  * GDPS( I,J ) * 100. / GRAV
*
               GDIE ( I,J,K )
     &                = CP * GDT( I,J,K )
     &                  * GDPS( I,J ) * 100. / GRAV
*
               GDPE ( I,J,K )
     &                = GDPE( I,J,K ) 
     &                  * GDPS( I,J ) * 100. / GRAV
*
               GDLE ( I,J,K )
     &                = EL * GDQ( I,J,K )
     &                  * GDPS( I,J ) * 100. / GRAV
*
               GDTE ( I,J,K )
     &                = GDKE ( I,J,K ) + GDIE ( I,J,K )
     &                + GDPE ( I,J,K ) + GDLE ( I,J,K )
*
               GDENS( I,J,K )
     &                = GDVOR ( I,J,K ) **2  
     &                  * GDPS( I,J ) * 100. / GRAV
*
 1100 CONTINUE
*
      DO 1200 I = 1, IMAX
         DO 1200 J = 1, JMAX
            GDTM ( I,J ) = GDPS ( I,J ) * 100. / GRAV
 1200 CONTINUE
*
      CALL HISTIN ( GDTM  , 'TM'  )
      CALL HISTIN ( GDKE  , 'KE'  )
      CALL HISTIN ( GDIE  , 'IE'  )
      CALL HISTIN ( GDPE  , 'PE'  )
      CALL HISTIN ( GDLE  , 'LE'  )
      CALL HISTIN ( GDTE  , 'TE'  )
      CALL HISTIN ( GDENS , 'ENS' )
*
      RETURN
      END
***********************************************************************
      SUBROUTINE DDIHRG      
*
*   [INTERNAL PARAM]
      INTEGER    ISTYPL
      DATA       ISTYPL / 1 /
*   [INTERNAL WORK]
      ReaL*8       VMISS
*
      CALL GZDBGT ( 'MISS' , VMISS )
*
      CALL     HISTRG
     I       ( 'TM    ', 'total mass          ' ,'kg    ', 'ASFC',
     I           VMISS , VMISS , VMISS , VMISS , ISTYPL  ,
     I         'CON'   , '   ' , 0     , 0     , 'XYZ'   ,'(1PE12.3)' )
*
      CALL     HISTRG
     I       ( 'KE    ', 'kinetic energy      ' ,'J     ', 'ALEV',
     I           VMISS , VMISS , VMISS , VMISS , ISTYPL  ,
     I         'CON'   , '   ' , 0     , 0     , 'XYZ'   ,'(1PE12.3)' )
*
      CALL     HISTRG
     I       ( 'IE    ', 'internal energy      ' ,'J     ', 'ALEV',
     I           VMISS , VMISS , VMISS , VMISS , ISTYPL  ,
     I         'CON'   , '   ' , 0     , 0     , 'XYZ'   ,'(1PE12.3)' )
*
      CALL     HISTRG
     I       ( 'PE    ', 'potential energy     ' ,'J     ', 'ALEV',
     I           VMISS , VMISS , VMISS , VMISS , ISTYPL  ,
     I         'CON'   , '   ' , 0     , 0     , 'XYZ'   ,'(1PE12.3)' )
*
      CALL     HISTRG
     I       ( 'LE    ', 'latent energy        ' ,'J     ', 'ALEV',
     I           VMISS , VMISS , VMISS , VMISS , ISTYPL  ,
     I         'CON'   , '   ' , 0     , 0     , 'XYZ'   ,'(1PE12.3)' )
*
      CALL     HISTRG
     I       ( 'TE    ', 'total energy         ' ,'J     ', 'ALEV',
     I           VMISS , VMISS , VMISS , VMISS , ISTYPL  ,
     I         'CON'   , '   ' , 0     , 0     , 'XYZ'   ,'(1PE12.3)' )
*
      CALL     HISTRG
     I       ( 'ENS   ', 'enstrophy            ' ,'kg    ', 'ALEV',
     I           VMISS , VMISS , VMISS , VMISS , ISTYPL  ,
     I         'CON'   , '   ' , 0     , 0     , 'XYZ'   ,'(1PE12.3)' )
*
      RETURN
      END
***********************************************************************
      SUBROUTINE SPCTLD      
     I         ( WDATA , HTTL  , KDIMD , NMO   )
*




* PACKAGE  ZXT10L32 !" 格子点数，波数 X-Z
*
      INTEGER     IMAX           !" 経度格子数
      INTEGER     JMAX           !" 緯度格子数
      INTEGER     IDIM           !" 経度配列サイズ
      INTEGER     JDIM           !" 緯度配列サイズ
      INTEGER     MMAX           !" 最大東西波数
      INTEGER     LMAX           !" 最大南北波数
      INTEGER     NMAX           !" 最大全波数
      INTEGER     MINT           !" 東西波数間隔
      INTEGER     KMAX           !" 鉛直層数
      INTEGER     KDIM           !" 鉛直配列サイズ
      INTEGER     NMDIM          !" 水平波数成分数
      INTEGER     JMXHF          !" JMAX/2+1
      INTEGER     NWORK          !" ワーク領域大きさ
*
      PARAMETER ( IMAX  = 32,     JMAX =  1,
     &            IDIM  = IMAX+1, JDIM = JMAX,
     &            MMAX  = 10,     LMAX =  0,    NMAX = MMAX, MINT = 1,
     &            KMAX  = 32,     KDIM = KMAX,
     &            NMDIM = (MMAX/MINT+1)*(2*(NMAX+1)-MMAX)
     &                  - (NMAX-LMAX)/MINT*(NMAX-LMAX+1),
     &            JMXHF = JMAX/2+1,
     &            NWORK = IDIM*JDIM*KMAX*12                       )

*
*   [INPUT]
      INTEGER    KDIMD
      ReaL*8       WDATA ( NMDIM, KDIMD )
      CHARACTER  HTTL  *(*)
      INTEGER    NMO   ( 2, 0:MMAX, 0:LMAX ) 
*
*   [INTERNAL WORK]
      ReaL*8       WAMP  ( 0:NMAX, KMAX+1 )
      LOGICAL    ODIAG
      INTEGER    JFLDIA, NCOLS, ISTR, IEND, IINT
      INTEGER    K,L,M,N
*
      CALL QDIAGP
     O     ( ODIAG , JFLDIA, NCOLS,
     O       ISTR  , IEND  , IINT ,
     I       KDIMD                  )
*
      IF ( .NOT. ODIAG )       RETURN
*
      CALL RESET ( WAMP, (NMAX+1)*KDIMD )
*
      DO 2100 K = 1, KDIMD
         DO 2110 N = 0 , NMAX
            DO 2120 M = 0 , MMAX, MINT
               L = N - M
               IF ( ( L .GE. 0 ).AND.( L .LE. LMAX ) ) THEN
                  WAMP ( N,K ) = WAMP ( N,K )
     &                         + WDATA( NMO(1,M,L),K ) ** 2
     &                         + WDATA( NMO(2,M,L),K ) ** 2
               ENDIF
 2120      CONTINUE
 2110   CONTINUE
 2100 CONTINUE
*
      WRITE ( JFLDIA, * ) '@@@@@@@@@@@@< SPECTRAL AMPLITUDE ', HTTL,
     &                    '>@@@@@@@@@@@@@'
*
      DO 3100 K = 1, KDIMD
         CALL PRINTD
     I       ( JFLDIA, '(1PE12.4)', WAMP  ,
     I         1     , NMAX+1 , 1     , K    ,
     I         NMAX+1, KDIMD  , NCOLS         )
 3100 CONTINUE
*
      RETURN
      END
