* PACKAGE DDIAG  !" 力学 診断 (弱非線形解析)
*
*"  [HIS]  90/05/19(numaguti)
*"         96/09/14(takepiro)  for baro, weak non-linear analysis
*
***********************************************************************
      SUBROUTINE DDIAG      !"  診断出力 
     I         ( GDU   , GDV   , GDVOR  )
*
#ifdef SYS_IBMS
      INCLUDE   (ZCDIM)                      !" 格子点数, 波数
      INCLUDE   (ZCCOM)                      !" 文字列文字数
#else
#include        "zcdim.F"                    !" 格子点数, 波数
#include        "zccom.F"                    !" 文字列文字数
#endif
*
*   [INPUT]
      REAL       GDU   ( IDIM, JDIM, KMAX )  !" 西風  Ｕ
      REAL       GDV   ( IDIM, JDIM, KMAX )  !" 南風  Ｖ
      REAL       GDVOR ( IDIM, JDIM, KMAX )  !" 渦度  ζ
*
*   [INTERNAL WORK]
      COMMON    /COMWRK /
     &           GDKE  , GDPV , GDLENS , GDAM
*
      REAL       GDKE  ( IDIM, JDIM, KMAX )  !" 運動エネルギー
      REAL       GDAM  ( IDIM, JDIM, KMAX )  !" 角運動量
      REAL       GDPV  ( IDIM, JDIM, KMAX )  !" ポテンシャル渦度
      REAL       GDENS ( IDIM, JDIM, KMAX )  !" エンストロフィー
*
      REAL       ALON  ( IDIM )              !" 経度
      REAL       DLON  ( IDIM )              !" 経度荷重
      REAL       ALAT  ( JDIM )              !" 緯度
      REAL       DLAT  ( JDIM )              !" 緯度荷重
*
      REAL       CORIOL ( IDIM, JDIM )       !" コリオリ係数 ｆ
*
      REAL       GDWA  ( IDIM, JDIM, KMAX )  !" wave activity
      REAL       GDBETA( IDIM, JDIM, KMAX )  !" 帯状平均絶対渦度の緯度微分
*
      REAL       GDVORD ( IDIM, JDIM, KMAX ) !" 擾乱相対渦度  ζ'
      REAL       GDVORM ( JDIM, KMAX )       !" 帯状平均相対渦度 \bar{ζ}
      REAL       GDTVM  ( JDIM, KMAX )       !" 帯状平均絶対渦度 \bar{ζ+f}
      REAL       WDBETA ( NMDIM , KMAX )     !" 帯状平均絶対渦度 \bar{ζ+f}
*
      REAL       VMISS
      INTEGER    I,J,K
      LOGICAL    OFIRST
      DATA       OFIRST / .TRUE. /
*
      SAVE       OFIRST, ALAT, DLAT, CORIOL, VMISS
*
      IF ( OFIRST ) THEN
         OFIRST = .FALSE.
         CALL ACLON
     O           ( ALON , DLON )
*
         CALL ACLAT
     O           ( ALAT , DLAT )
*
         CALL DCCORL
     O        ( CORIOL, 
     I          ALAT  , DLAT )
*         
         CALL DHSTRG
         CALL GZDBGT( 'MISS', VMISS )
      ENDIF
*
*"         < 1. エネルギー，エンストロフィー >
*
      DO 1100 K = 1, KMAX
         DO 1100 J = 1, JDIM
            DO 1100 I = 1, IDIM
*
               GDKE ( I,J,K )
     &                = ( GDU( I,J,K )**2 + GDV( I,J,K )**2 )/ 2.  
*
               GDAM ( I,J,K )
     &                = GDU( I,J,K ) * ER * COS( ALAT(J) )
*
               GDPV ( I,J,K )
     &                = GDVOR( I,J,K ) + CORIOL( I, J )
*
               GDENS( I,J,K )
     &                = GDVOR ( I,J,K ) **2 / 2.
*
 1100 CONTINUE
*
      CALL HISTIN ( GDKE  , 'KE'  )
      CALL HISTIN ( GDAM  , 'AM'  )
      CALL HISTIN ( GDPV  , 'PV'  )
      CALL HISTIN ( GDENS , 'ENS' )
*
*
*"          < 2. 弱非線形解析  > 
*
      CALL RESET( GDVORM, JDIM*KMAX )
      CALL RESET( GDTVM , JDIM*KMAX )
*
      DO 2100 K = 1, KMAX
         DO 2100 J = 1, JDIM
            DO 2110 I = 1, IDIM
               GDVORM( J, K ) = GDVORM( J, K ) 
     &              + GDVOR( I,J,K )*DLON(I)            !" 帯状平均相対渦度
               GDTVM( J, K ) = GDTVM( J, K )
     &              + ( CORIOL( I,J ) + GDVOR( I,J,K ) ) 
     &              * DLON(I)                           !" 帯状平均絶対渦度
 2110       CONTINUE 
*
            DO 2120 I = 1, IDIM
               GDVORD( I, J, K ) 
     &              = GDVOR( I,J,K ) - GDVORM( J,K )    !" 擾乱相対渦度
               GDBETA( I, J, K ) = GDTVM( J,K )
 2120       CONTINUE 
 2100 CONTINUE
*
      CALL G2W                  !" 平均絶対渦度(スペクトル)
     O         ( WDBETA,
     I           GDBETA,
     I           ' ', 'POS ',  KMAX  )
*
      CALL W2G                  !" 平均絶対渦度の y 微分
     O         ( GDBETA,
     I           WDBETA,
     I           'YGRA', 'POS ', KMAX )
*
      DO 2200 K = 1, KMAX
         DO 2200 J = 1, JDIM
            DO 2200 I = 1, IDIM
               IF ( GDBETA( I,J,K ) .NE. 0.0 ) THEN
                  GDWA( I, J, K ) 
     &                 = GDVORD( I,J,K )**2 * COS( ALAT(J) )**4
     &                 / ( 2 * GDBETA( I,J,K ) ) * ER**2
               ELSE
                  GDWA( I, J, K ) = VMISS
               ENDIF
 2200 CONTINUE 
*
      CALL HISTIN ( GDWA  , 'WA'  )
*
      RETURN
      END
***********************************************************************
*"          << 標準時間平均出力の登録 >>
***********************************************************************
      SUBROUTINE DHSTRG
*
      REAL       VMISS
      INTEGER    ISTYPL
      DATA       ISTYPL / 1 /
*
*"         < 3. ヒストリー登録 >
*
      CALL GZDBGT( 'MISS', VMISS )
*
      CALL     HISTRG
     I       ( 'WA     ', 'wave activity       ' ,'m2/s  ', 'ASFC',
     I           VMISS , VMISS , VMISS , VMISS , ISTYPL  ,
     I         'CON'   , '   ' , 0     , 0     , 'X'     ,'(F12.3)' )
*
      RETURN
      END
