* PACKAGE A1DIM      !" 鉛直1次元メイン
*
*"  [HIS] 90/08/31(numaguti)
*         93/09/24(ishiwatari)
*
***********************************************************************
      PROGRAM GCM1D     !" 鉛直一次元テストモデル
*
*
*"       << 鉛直一次元モデル >>
*
*"                 Version 0.50              90/08/24 ( 沼口 敦 )
*
*   [PARAM] 
#if   HITAC
      INCLUDE   (ZCDIM)
      INCLUDE   (ZHDIM)
      INCLUDE   (ZCCOM)                      !" 標準物理定数
#else
#include        "zcdim.F"
#include        "zhdim.F"
#include        "zccom.F"                    !" 標準物理定数
#endif
*
*"  [VAR]
      REAL       GBU   ( IDIM, JDIM, KMAX )  !" 西風  ｕ
      REAL       GBV   ( IDIM, JDIM, KMAX )  !" 南風  ｖ
      REAL       GBT   ( IDIM, JDIM, KMAX )  !" 温度  Ｔ
      REAL       GBPS  ( IDIM, JDIM       )  !" 地表気圧
      REAL       GBQ   ( IDIM, JDIM, KMAX )  !" 比湿  ｑ
      REAL       GBVOR ( IDIM, JDIM, KMAX )  !" 渦度  
      REAL       GBDIV ( IDIM, JDIM, KMAX )  !" 発散  Ｄ
      REAL       GBTG  ( IDIM, JDIM )        !" 地面温度Ｔｇ
*"       : 格子点データ
*
      REAL       GTUA  ( IDIM, JDIM, KMAX )  !" 東西運動量変化項ＵＡ
      REAL       GTVA  ( IDIM, JDIM, KMAX )  !" 南北運動量変化項ＶＡ
      REAL       GTH   ( IDIM, JDIM, KMAX )  !" 温度時間変化項  Ｈ
      REAL       GTR   ( IDIM, JDIM, KMAX )  !" 比湿時間変化項  Ｒ
      REAL       GTTG  ( IDIM, JDIM )        !" 地表温度変化率
*"       : 格子点時間変化項
*
      INTEGER    ISTEP                       !" 通しステップ数
      INTEGER    ITA                         !" 通し時間(t)
      INTEGER    ITB                         !" 通し時間(t-Δt)
      INTEGER    IDATEA( 6 )                 !" 時刻：年月日時分秒
      INTEGER    IDATEB( 6 )                 !" 時刻：年月日時分秒
      REAL       DELT                        !" 時間刻みΔt
      REAL       DELTP                       !" 物理過程時間刻み 2 Δt
*                                               !" 石渡 93/09/24 追加
      LOGICAL    OADVNC                      !" 時刻が進行するか否か
      LOGICAL    ORSTR                       !" 再出発出力するか否か
*"        : 時刻変数
*
      DATA       ISTEP        / 0  /
      DATA       ITA, ITB     / 0 , 0 /
*
*   [INTERNAL PARM] 
      INTEGER    IDATEI ( 6 )
      DATA       IDATEI       / 1992, 1, 1, 0, 0, 0 /
      NAMELIST   / NMDATE / IDATEI 
*
      INTEGER    ITSTRT                      !" 計算開始時刻
      INTEGER    ITEND                       !" 計算終了時刻
      INTEGER    IORSTR                      !" 出力間隔：再出発
*
      INTEGER    IDELT                       !" 標準時間刻み
      INTEGER    INISTP                      !" 初期ステップ数
      REAL       TFIL                        !" 時間フィルター係数
*
*"    COMMON /COMCON/ ( include ZCCOM )
*     REAL       CP                          !" 大気定圧比熱
*     REAL       EL                          !" 水の凝結の潜熱
*     REAL       ER                          !" 地球半径
*     REAL       GRAV                        !" 重力加速度
*     REAL       RAIR                        !" 大気気体定数
*     REAL       RVAP                        !" 水蒸気気体定数
*     REAL       CPVAP                       !" 水蒸気定圧比熱
*     REAL       DH2O                        !" 水の密度
*     REAL       EPSV                        !" 水と大気の分子量比
*     REAL       ES0                         !" ０℃飽和蒸気圧
*     REAL       STB                         !" Stefan-Boltzman
*"    COMMON /COMCON/ end
*
*   [INTERNAL WORK] 
      REAL       ALAT  ( JDIM )              !" 緯度
      REAL       DLAT  ( JDIM )              !" 緯度荷重
      REAL       ALON  ( IDIM )              !" 経度
      REAL       DLON  ( IDIM )              !" 経度荷重
*
      REAL       SIG   ( KMAX )              !" σレベル(整数)
      REAL       SIGM  (KMAX+1)              !" σレベル(半整数)
      REAL       DSIG  ( KMAX )              !" Δσ(整数)
      REAL       DSIGM (KMAX+1)              !" Δσ(半整数)
*
      REAL       CORIOL( JDIM )              !" コリオリ係数 ｆ
      REAL       GGU   ( IDIM, JDIM, KMAX )  !" 地衡風 ｕg
*
      REAL       PI, OMEGA
      INTEGER    IFPAR, JFPAR, IERR
      INTEGER    I, J, K, II
*
*   [WORK] 
      COMMON    /COMWRK/ WORK
      REAL       WORK  ( NWORK )             !" ワーク領域
*
      CALL YPREP
*
*"           <<<<<<<<<<<<< 初期設定 >>>>>>>>>>>>>>
*
*
*"         < 1. 実験パラメーターのセット >
*
      OPEN     ( 1, FILE='PARM1D', IOSTAT=IERR,
     &           FORM='FORMATTED',
     &           ACCESS='SEQUENTIAL'                 )
      CALL       SETNML ( 1, 6 )
*
      CALL       SETPAR
     O         ( ITSTRT, ITEND , IORSTR,
     O           IDELT , INISTP, TFIL            )
*
      CALL PCONST
*
*"         < 2. モデル座標定数の計算 >
*
      CALL SETCOR
     O         ( ALON  , DLON  ,
     O           ALAT  , DLAT  ,
     O           SIG   , DSIG  ,
     O           SIGM  , DSIGM  )
*
      PI    = ATAN2 ( 0.  , -1. )
      OMEGA = 2.   * PI / 86400.
      DO 200 J = 1, JDIM
         CORIOL ( J ) = 2. * OMEGA * SIN( ALAT(J) )
  200 CONTINUE
*
*"         < 5. 初期値の生成 >
*
      CALL   REWNML ( IFPAR , JFPAR )
      READ   ( IFPAR, NMDATE, END=480 )
 480  WRITE  ( JFPAR, NMDATE )
      DO 490 II = 1, 6
         IDATEA( II ) = IDATEI( II )
         IDATEB( II ) = IDATEI( II )
 490  CONTINUE
*
      DO 500 K = 1, KMAX
         DO 500 J = 1, JDIM
            DO 500 I = 1, IDIM
               GGU  ( I,J,K ) = 20.  *(1.5-SIG(K))
               GBU  ( I,J,K ) = GGU ( I,J,K )
               GBV  ( I,J,K ) = 0.
               GBVOR( I,J,K ) = 0.
               GBDIV( I,J,K ) = 0.
               GBT  ( I,J,K ) = 200. + SIG(K) * 100.
               GBQ  ( I,J,K ) = 1.E-5
  500 CONTINUE
      DO 510 J = 1, JDIM
         DO 510 I = 1, IDIM
               GBPS( I,J ) = 1000.  
               GBTG( I,J ) = GBT ( I,J,1 )
  510 CONTINUE
*
      WRITE ( 6,* ) ' ###### START ######## ', IMAX, JMAX, KMAX
 5000 CONTINUE
*"                                               <===== repeat MAIN ===
*
*"         <<<<<<<<< メイン・ループ始まり >>>>>>>>
*
         ISTEP = ISTEP + 1
*
*"         < 6 時間制御 >
*
         CALL TIMSTP
     O         ( DELT  , OADVNC, ORSTR ,
     I           ITA   , ITB   , IDATEA, IDATEB, ISTEP  ,
     C           IDELT , INISTP, ITEND , IORSTR          )
*
         CALL SAVEGQ
     I         ( GBQ   )
*
*"         < 7. 物理過程 >
*"                           93/09/24 石渡 変更
*
         DELTP = 2 * DELT
*
         CALL PHYSCS
     O         ( GTUA  , GTVA  , GTH   , GTR   ,
     M           GBTG  ,
     I           GBU   , GBV   , GBT   , GBPS  , GBQ   ,
     I           GBVOR , GBDIV ,
     I           ITB   , IDATEB, DELTP , DELT  ,
     C           ALON  , ALAT  ,
     C           SIG   , SIGM  , DSIG  , DSIGM          )
*
*"                < 6.2 地面温度の積分 >
*
         CALL TGINTG
     M         ( GBTG  ,
     I           GTTG  , DELT  )
*
*"         < 8. 力学過程 >
*
         DO 800 K = 1, KMAX
            DO 800 J = 1, JDIM
               DO 800 I = 1, IDIM
                  GBU  ( I,J,K ) = GBU ( I,J,K )
     &                           + GTUA( I,J,K ) * DELT
     &                           + GBV ( I,J,K )
     &                              * CORIOL ( J ) * DELT
*
                  GBV  ( I,J,K ) = GBV ( I,J,K )
     &                           + GTVA( I,J,K ) * DELT
     &                           + ( GGU ( I,J,K ) - GBU  ( I,J,K ) )
     &                              * CORIOL ( J ) * DELT
*
                  GBT  ( I,J,K ) = GBT ( I,J,K )
     &                           + GTH ( I,J,K ) * DELT
*
                  GBQ  ( I,J,K ) = GBQ ( I,J,K )
     &                           + GTR ( I,J,K ) * DELT
  800    CONTINUE
*
         CALL PADJST
     M         ( GBU   , GBV   , GBT   , GBPS  , GBQ   ,
     I           GBVOR , GBDIV ,
     I           DELT  ,
     C           SIG   , SIGM  , DSIG  , DSIGM          )
*
         CALL   DIAGP ( GBU   , 'U'   , KMAX  , '(F12.3)'  )
         CALL   DIAGP ( GBV   , 'V'   , KMAX  , '(F12.3)'  )
         CALL   DIAGP ( GBT   , 'T'   , KMAX  , '(F12.3)'  )
         CALL   DIAGP ( GBQ   , 'Q'   , KMAX  , '(3PF12.3)')
         CALL   DIAGP ( GBTG  , 'TG'  , 1     , '(F12.3)'  )
*
         CALL    ADVSTP
     O         ( ITA   , ITB   , IDATEA, IDATEB,
     I           OADVNC                          )
*
      IF ( ITB .LT. ITEND ) GOTO 5000
*
*"          <<<<<<<< メイン・ループ終わり >>>>>>>>
*
*"                                               =====> until MAIN ===>
      CALL YFINE
*
      STOP
      END
************************************************************************
      SUBROUTINE RDCON
     O         ( GDATA , ITNEXT,
     I           HFILE , HITEM , DEFVAL , HCORD , ITSEL  )
*
*   [PARAM] 
#if   HITAC
      INCLUDE   (ZCDIM)                      !" 格子点数，波数
#elif SX3
      INCLUDE   "zcdim.F"                    !" 格子点数，波数
#else
#include        "zcdim.F"                    !" 格子点数，波数
#endif
*
*   [OUTPUT] 
      REAL       GDATA ( * )                 !" データ
      INTEGER    ITNEXT
*
*   [INPUT] 
      CHARACTER  HFILE *(*)                  !" ファイル名
      CHARACTER  HITEM *(*)                  !" 識別名称
      REAL       DEFVAL                      !" デフォルト値
      CHARACTER  HCORD *(*)                  !" 座標
      INTEGER    ITSEL
*
*   [INTERNAL PARM] 
      INTEGER    ILARGE
      DATA       ILARGE / 2147483647 /       !" 整数最大値
*
*   [INTERNAL WORK] 
      INTEGER    KDIMD, NDIMD
*
      IF      ( HCORD .EQ. 'ALEV' ) THEN
         KDIMD = KMAX
      ELSE IF ( HCORD .EQ. 'AMLEV' ) THEN
         KDIMD = KMAX+1
      ELSE
         KDIMD = 1
      ENDIF
      NDIMD = IDIM * JDIM * KDIMD
*
      CALL  SETV ( GDATA, DEFVAL, NDIMD  )
*
      ITNEXT = ILARGE
*
      RETURN
      END
************************************************************************
      SUBROUTINE HISTIN
     I         ( GDATA ,
     C           HITEM  )
*
*"          << ダミー >>
*
*   [INPUT] 
      REAL       GDATA
      CHARACTER  HITEM *(*)
*
      RETURN
      END
************************************************************************
      SUBROUTINE HISTRG
     I         ( HITEM , HTITL , HUNIT , HCORD ,
     I           DMINX , DMAXX , DIVSX , DIVLX , ISTYPX,
     I           FILEM , TUNITM, INTOUM, INTSMM, AVRTYM, DFORMM )
*
*"          << ダミー >>
*
      CHARACTER  HITEM *(*), HTITL *(*), HUNIT *(*), HCORD *(*)
      REAL       DMINX, DMAXX, DIVSX, DIVLX
      INTEGER    ISTYPX
      CHARACTER  FILEM *(*), TUNITM*(*)
      INTEGER    INTOUM, INTSMM
      CHARACTER  AVRTYM*(*), DFORMM*(*)
*
      RETURN
      END
************************************************************************
      SUBROUTINE SAVEGQ
     I         ( GBQ    )
*
#if   HITAC
      INCLUDE   (ZCDIM)                      !" 格子点数, 波数
#else
#include        "zcdim.F"                    !" 格子点数, 波数
#endif
*
*   [INPUT] 
      REAL       GBQ   ( IDIM*JDIM*KMAX )    !" 比湿  ｑ(t-Δt)
*
*   [INTERNAL SAVE] 
      REAL       GBBQ  ( IDIM*JDIM*KMAX )    !" 比湿  ｑ セーブ値
      SAVE
*
      CALL COPY ( GBBQ  , GBQ   , IDIM*JDIM*KMAX )
*
      RETURN
*======================================================================
      ENTRY      LOADGQ
     O         ( GBQ    )
*
      CALL COPY ( GBQ   , GBBQ  , IDIM*JDIM*KMAX )
*
      RETURN
      END
***********************************************************************
      SUBROUTINE GAUSS
     O         ( QSINLA , QGW    ,
     I           JMAX             )
*
*"          << ダミー >>
*
      INTEGER    JMAX
      REAL       QSINLA
      REAL       QGW

      RETURN
      END
***********************************************************************
      SUBROUTINE IFLOPN
     O         ( IFILE , IERR  ,
     I           HFILE , IFILED, HACT  , HFORM   )
*
*   [OUTPUT] 
      INTEGER    IFILE
      INTEGER    IERR
*
*   [INPUT] 
      CHARACTER  HFILE *(*)
      INTEGER    IFILED
      CHARACTER  HACT  *(*)
      CHARACTER  HFORM *(*)
*
      IFILE = 6
      IERR  = 0
*
      RETURN
      END
