* PACKAGE P2ADMN  !" 物理過程(2) メインパート
*
*"履歴 92/04/27 沼口 敦             重力波抵抗(masaaki)の導入
*"履歴 93/08/29 保坂征宏            長波放射変更
*"履歴 94/04/28 石渡正樹            平衡モデル用変更
*"
***********************************************************************
      SUBROUTINE PHYSCS      !" 物理過程メイン
     O         ( GTU   , GTV   , GTT   , GTQ   , 
     M           GDTG  ,
     I           GDU   , GDV   , GDT  , GDPS  , GDQ   ,
     I           GDVOR , GDDIV , 
     I           IT    , IDATE , DELT  , DELTM , 
     C           ALON  , ALAT  ,
     C           SIG   , SIGM  , DSIG  , DSIGM          )
*
*"  物理過程： 多色放射  鉛直拡散 地表 
*
*   [PARAM] 
#ifdef SYS_IBMS
      INCLUDE   (ZCDIM)                      !" 格子点数，波数
      INCLUDE   (ZCCOM)                      !" 標準物理定数
#else
#include        "zcdim.F"                    !" 格子点数，波数
#include        "zccom.F"                    !" 標準物理定数
#endif
*
*   [OUTPUT] 
      REAL       GTU   ( IDIM, JDIM, KMAX )  !" ＵＡ(ダミー)
      REAL       GTV   ( IDIM, JDIM, KMAX )  !" ＶＡ(ダミー)
      REAL       GTT   ( IDIM, JDIM, KMAX )  !" 温度時間変化項 Ｈ
      REAL       GTQ   ( IDIM, JDIM, KMAX )  !" 比湿時間変化項 Ｒ
      REAL       GTTG  ( IDIM, JDIM )        !" 地表温度変化率
*
*   [INPUT] 
      REAL       GDU   ( IDIM, JDIM, KMAX )  !" 西風  ｕ(ダミー)
      REAL       GDV   ( IDIM, JDIM, KMAX )  !" 南風  ｖ(ダミー)
      REAL       GDT   ( IDIM, JDIM, KMAX )  !" 温度  Ｔ
      REAL       GDPS  ( IDIM, JDIM       )  !" 地表気圧 Ｐｓ
      REAL       GDQ   ( IDIM, JDIM, KMAX )  !" 比湿  ｑ
      REAL       GDVOR ( IDIM, JDIM, KMAX )  !" 渦度  ζ(ダミー)
      REAL       GDDIV ( IDIM, JDIM, KMAX )  !" 発散  Ｄ(ダミー)
      REAL       GDTG  ( IDIM, JDIM )        !" 地表面温度
*
      INTEGER    IT                          !" 通し時間
      INTEGER    IDATE ( * )                 !" 時刻：年月日時分秒
      REAL       DELT                        !" 物理過程時間刻み 2 Δt
      REAL       DELTM                       !" モデル  時間刻み   Δt
*                                                !" 保坂 93/08/29 追加
      REAL       ALON  ( IDIM )              !" 経度
      REAL       ALAT  ( JDIM )              !" 緯度
      REAL       SIG   ( KMAX )              !" σレベル(整数)
      REAL       SIGM  (KMAX+1)              !" σレベル(半整数)
      REAL       DSIG  ( KMAX )              !" Δσ
      REAL       DSIGM (KMAX+1)              !" Δσ(半整数)
*
*   [INTERNAL SAVE] 
      REAL       RFLUX  ( IDIM, JDIM, KMAX+1 )   !" 長波フラックス
      REAL       RFLUXC ( IDIM, JDIM, KMAX+1 )   !" 補正した長波フラックス
*                                                   !" (ダミー)
*
      REAL       SFLUX  ( IDIM, JDIM, KMAX+1 )   !" 短波フラックス
      REAL       DRFLUX ( IDIM, JDIM, KMAX+1, 0:1 ) !" 長波温度変化
*
      REAL       GRHC   ( IDIM, JDIM )        !" 地表熱容量
      REAL       GRALB  ( IDIM, JDIM )        !" 地表アルベド
      REAL       GRBET  ( IDIM, JDIM )        !" 地表湿潤度
      REAL       GRZ0   ( IDIM, JDIM )        !" 地表粗度長
      REAL       GRZSD  ( IDIM, JDIM )        !" 地表高度分散
      INTEGER    IDSFC  ( IDIM, JDIM )        !" 地表状態
      REAL       GFLUX  ( IDIM, JDIM )        !" 地中熱フラックス
*
      REAL       UFLUX  ( IDIM, JDIM, KMAX+1)      !" Ｕのフラックス
      REAL       VFLUX  ( IDIM, JDIM, KMAX+1)      !" Ｖのフラックス
      REAL       TFLUX  ( IDIM, JDIM, KMAX+1)      !" Ｔのフラックス
      REAL       QFLUX  ( IDIM, JDIM, KMAX+1)      !" ｑのフラックス
      SAVE
*
*   [INTERNAL PARM] 
      REAL       P0                       !" 相対湿度
      DATA       P0 / 1000.0 /
*
      NAMELIST  /NMADEQ/ P0
*
      LOGICAL    OPBDGT                       !" 収支を出力
      LOGICAL    OFIRST
      DATA       OPBDGT / .TRUE. /
      DATA       OFIRST / .TRUE. /
      SAVE OFIRST, OPBDGT
*
*
*   [INTERNAL WORK] 
      REAL       GDTM   ( IDIM, JDIM, KMAX+1)  !" 温度：半整数
      REAL       GDP    ( IDIM, JDIM, KMAX )   !" 気圧  Ｐ
      REAL       GDPM   ( IDIM, JDIM, KMAX+1 ) !" 気圧（半整数）
      REAL       GDZ    ( IDIM, JDIM, KMAX )   !" 高度
      REAL       GDZM   ( IDIM, JDIM, KMAX+1 ) !" 高度（半整数）
      REAL       GTTRL  ( IDIM, JDIM, KMAX )   !" 温度時間変化：長波
      REAL       GTTRS  ( IDIM, JDIM, KMAX )   !" 温度時間変化：短波
      REAL       GDQG   ( IDIM, JDIM )   !" 地表q
*
      REAL       GQPR   ( IDIM, JDIM       )  !" 降水量
      REAL       GQCND  ( IDIM, JDIM, KMAX )  !" 凝結加熱率
      REAL       GSCND  ( IDIM, JDIM, KMAX )  !" 凝結比湿変化
      REAL       GSRMNQ ( IDIM, JDIM, KMAX )  !" 負の比湿除去
      REAL       TRPPM  (  IDIM, JDIM )       !" 圏界面圧力
      REAL       RUFLUX  ( IDIM, JDIM, KMAX+1 )   !" 長波フラックス
      REAL       RDFLUX  ( IDIM, JDIM, KMAX+1 )   !" 長波フラックス
      INTEGER    TRPLEV ( IDIM, JDIM )            !" 圏界面の高さ
*
      REAL       RMTXS  ( IDIM, JDIM,      -1:1 ) !" Ｔ陰解行列(ダミー)
*
      INTEGER    I, J, K
      INTEGER    IFPAR, JFPAR
*      
*   [INTERNAL FUNC]
#ifdef SYS_IBMS
      INCLUDE   (ZQSAT)                      !" 飽和水蒸気圧
#else
#include        "zqsat.F"                    !" 飽和水蒸気圧
#endif
*
*   [ONCE] 
      IF ( OFIRST ) THEN
         WRITE ( 6,* ) ' PHYSICS PACKAGE EQUIB DATE=94/04/29'
         OFIRST = .FALSE.
         CALL PHSTRG
*
         CALL   REWNML ( IFPAR , JFPAR )
         READ   ( IFPAR, NMADEQ, END=190 )
  190    WRITE  ( JFPAR, NMADEQ )
*
      ENDIF
*
*"          < -1. 初期値設定 >
*
*"         こんなところで初期値設定 Part I
*"         なんとかせい
*
      GDT  ( 1,1,1 ) = 2.86207E+02
      GDTG ( 1,1   ) = 2.86207E+02
      GDT  ( 2,1,1 ) = 3.00196E+02
      GDTG ( 2,1   ) = 3.00196E+02
      GDT  ( 3,1,1 ) = 3.13129E+02
      GDTG ( 3,1   ) = 3.13129E+02
      GDT  ( 4,1,1 ) = 3.17744E+02
      GDTG ( 4,1   ) = 3.17744E+02
      GDT  ( 5,1,1 ) = 3.20652E+02
      GDTG ( 5,1   ) = 3.20652E+02
      GDT  ( 6,1,1 ) = 3.48713E+02
      GDTG ( 6,1   ) = 3.48713E+02
      GDT  ( 7,1,1 ) = 3.52171E+02
      GDTG ( 7,1   ) = 3.52171E+02
      GDT  ( 8,1,1 ) = 3.63477E+02
      GDTG ( 8,1   ) = 3.63477E+02
      DO 500 I=9,IDIM
         DO 500 J=1,JDIM
            GDT  ( I,J,1 ) = 315.361
            GDTG ( I,J   ) = 315.361
 500  CONTINUE 

      DO 1000 I=1,IDIM
         DO 1000 J=1,JDIM
c$$$            GDT  ( I,J,1 ) = 250. + REAL(I-1) * 10.0
c$$$            GDTG ( I,J ) = 250. + REAL(I-1) * 10.0
            GDPS ( I,J   ) = P0
            GDQG ( I,J   ) = FQSAT( GDTG (I,J ), GDPS( I,J )*100. )
            GDPS ( I,J   ) = GDPS( I,J )*( 1.0 + GDQG( I,J )/EPSV )
c            GDPS ( I,J ) = GDPS( I,J )*( GDQG( I,J )/EPSV )
            GDQG ( I,J   ) = FQSAT( GDTG (I,J ), GDPS( I,J )*100. )
            GDQ  ( I,J,1 ) = FQSAT( GDTG (I,J ), GDPS( I,J )*100. )
c            GDQG( I,J ) = FQSAT( GDTG (I,J ), GDPS( I,J )*100. )
c            GDQ( I,J,1 ) = FQSAT( GDT (I,J,1 ), GDP( I,J,1 )*100. )
 1000 CONTINUE 
*
*"         < 0. 温度，気圧，高度 >
*
      CALL TINTP
     O         ( GDTM  ,
     I           GDT   ,
     C           SIG   , SIGM  )
*
      CALL PZLEVS
     O         ( GDP   , GDPM , GDZ   , GDZM  ,
     I           GDPS  , GDT  , GDTM  ,
     C           SIG   , SIGM , DSIG  , DSIGM   )
*
*"         < 1. 地表モデル >
*
      CALL PGRNDM
     M         ( GDTG  ,
     M           GRHC  , GRALB , GRBET , GRZ0  , GRZSD , 
     M           IDSFC , GFLUX ,
     I           UFLUX , VFLUX , TFLUX , QFLUX , RFLUX , SFLUX ,
     I           GDU   , GDV   , GDTM  , GDPS  , GDQ   ,
     I           IT    , IDATE , DELT  ,
     C           ALON  , ALAT                                   )
*
*"         こんなところで初期値設定 Part II
*"         なんとかせい
*
      DO 1100 I=1,IDIM
         DO 1100 J=1,JDIM
c$$$            GDT  ( I,J,1 ) = 250. + REAL(I-1) * 10.0
c$$$            GDTG ( I,J ) = 250. + REAL(I-1) * 10.0
c            GDT( I,J,1 ) = 315.361
c            GDTG( I,J ) = 315.361
c$$$            GDQ( I,J,1 ) = FQSAT( GDT (I,J,1 ), GDP( I,J,1) )
 1100 CONTINUE 

      GDT  ( 1,1,1 ) = 2.86207E+02
      GDTG ( 1,1   ) = 2.86207E+02
      GDT  ( 2,1,1 ) = 3.00196E+02
      GDTG ( 2,1   ) = 3.00196E+02
      GDT  ( 3,1,1 ) = 3.13129E+02
      GDTG ( 3,1   ) = 3.13129E+02
      GDT  ( 4,1,1 ) = 3.17744E+02
      GDTG ( 4,1   ) = 3.17744E+02
      GDT  ( 5,1,1 ) = 3.20652E+02
      GDTG ( 5,1   ) = 3.20652E+02
      GDT  ( 6,1,1 ) = 3.48713E+02
      GDTG ( 6,1   ) = 3.48713E+02
      GDT  ( 7,1,1 ) = 3.52171E+02
      GDTG ( 7,1   ) = 3.52171E+02
      GDT  ( 8,1,1 ) = 3.63477E+02
      GDTG ( 8,1   ) = 3.63477E+02
      DO 600 I=9,IDIM
         DO 600 J=1,JDIM
            GDT  ( I,J,1 ) = 315.361
            GDTG ( I,J   ) = 315.361
 600  CONTINUE 
      DO 1200 I=1,IDIM
         DO 1200 J=1,JDIM
            GDQ( I,J,1 ) = FQSAT( GDT (I,J,1 ), GDP( I,J,1) )
 1200 CONTINUE 

*
*"         < 2. 湿潤断熱線 >
*
      CALL RESET ( GQPR  , IDIM*JDIM      )
      CALL RESET ( GQCND , IDIM*JDIM*KMAX )
      CALL RESET ( GSCND , IDIM*JDIM*KMAX )
      CALL RESET ( GSRMNQ, IDIM*JDIM*KMAX )
*
      CALL ADIABT
     M         ( GDT   , GDQ   ,
     I           GDP   , GDPM  , GDZ   ,
     I           GDTG  , GDQG  )
*
*"         < 3. 乾燥対流調節 >
*
      CALL DADJST
     M         ( GDT   ,
     I           GDP   , GDPM  , DELT   )
*
*"         < 4. 負の水蒸気除去 >
*
      CALL RMNEGQ
     M         ( GDQ   , GSRMNQ,
     I           GDPM  ,
     I           DELT           )
*
*"         < 5. 放射 flux >
*
      CALL HISTIN ( GDT , 'ADT'  )
      CALL HISTIN ( GDQ , 'ADQ'  )
*
      CALL TINTP
     O         ( GDTM  ,
     I           GDT   ,
     C           SIG   , SIGM  )
*
      CALL RADFLX                                !" 保坂 93/08/29 変更
     O         ( RFLUX , RMTXS , DRFLUX,
     O           SFLUX ,
     I           GDTM  , GDTG  , GDQ   , GDPM  , 
     I           IT    , IDATE ,
     C           ALON  , ALAT  , GRALB          )
*
         CALL UDFLUX
     O         ( RUFLUX , RDFLUX )
*
         write (6,*) 'in p2adeq;ij,k,rflux'
         do 2 i=1,idim
            do 2 k=1,kmax+1
               write (6,*) i,k,rflux(i,1,k)
 2          continue 
      CALL HISTIN ( RFLUX  , 'ALRAD'  )
      CALL HISTIN ( RUFLUX , 'AUFLUX'  )
      CALL HISTIN ( RDFLUX , 'ADFLUX'  )
*
*
*
*"         < 6. 温度構造決定 >
*
      CALL TRPPAS
     M        ( GDT   , GDQ,
     O          TRPLEV,
     I          GDTG  ,
     I          GDPM  , 
     I          RFLUX ,
     C          SIG   , SIGM )
*
      DO 6000 I=1,IDIM
         DO 6000 J=1,JDIM
            TRPPM(I,J) = SIGM(TRPLEV(I,J))*GDPM(I,J,1)
 6000 CONTINUE 
*
*"         < 7. 放射フラックス再計算 >
*
      CALL TINTP
     O         ( GDTM  ,
     I           GDT   ,
     C           SIG   , SIGM  )
*
      CALL RADFLX                                !" 保坂 93/08/29 変更
     O         ( RFLUX , RMTXS , DRFLUX,
     O           SFLUX ,
     I           GDTM  , GDTG  , GDQ   , GDPM  , 
     I           IT    , IDATE ,
     C           ALON  , ALAT  , GRALB          )
*
         CALL UDFLUX
     O         ( RUFLUX , RDFLUX )
*
*"         < 8. 放射による温度変化率 >
*
      CALL RADTND
     O         ( GTTRL , GTTRS ,
     I           RFLUX , SFLUX , GDPM   )       !" 保坂 93/08/29 記号変更
*
      write (6,*) 'after radtnd'
      do 1 i=1,idim
         do 1 k=1,kmax
            write(6,*) i,k,gttrl(i,1,k)
 1    continue 
*
*"         < 9. 出力 >
*
      CALL HISTIN ( TRPPM , 'TRPPM'  )
      CALL HISTIN ( GTTRL , 'QRADL' )
      CALL HISTIN ( GTTRS , 'QRADS' )
      CALL HISTIN ( RFLUX , 'LRAD'  )           !" 石渡 94/05/10
      CALL HISTIN ( RUFLUX , 'RUFLUX'  )
      CALL HISTIN ( RDFLUX , 'RDFLUX'  )
      CALL HISTIN ( RFLUX(1,1,KMAX+1),'OLR' )  !" 保坂 93/08/29 記号変更
      CALL HISTIN ( RFLUX,            'SLR' )  !" 保坂 93/08/29 記号変更
*
      IF ( OPBDGT ) THEN                        
         CALL PDIAG
     I         ( TFLUX , QFLUX , RFLUXC, SFLUX ,  !" 保坂 93/08/29 記号変更
     I           GTT   , GTQ   , GTTG  ,
     I           GDT   , GDQ   , GDTG  , GDPM  ,
     C           GRHC                           )
      ENDIF
*
      RETURN
      END
***********************************************************************
      SUBROUTINE PHSTRG    !" 物理過程時間平均出力の登録
*
      INTEGER    ISTYPL
      DATA       ISTYPL / 1 /
      REAL       VMISS
*
      CALL GZDBGT( 'MISS', VMISS )
*
      CALL HISTRG
     I       ( 'ADT'  , 'Temperature      ' ,'K    ', 'ALEV',
     I          VMISS     , VMISS , 2.5  , 10. , ISTYPL  ,
     I         'CON'   , '   ' , 0     , 0     , 'X'     ,'(F12.3)' )
*
      CALL     HISTRG
     I       ( 'ADQ   ', 'specific humidity   ' ,'g/g   ', 'ALEV',
     I           VMISS , VMISS , 1.E-3 , 5.E-3 , ISTYPL  ,
     I         'CON'   , '   ' , 0     , 0     , 'X'     ,'(3PF12.3)' )
*
      CALL HISTRG
     I       ( 'RAIN'  , 'precipitation      ' ,'W/m**2', 'ASFC',
     I          0.     , VMISS , 200.  , 1000. , ISTYPL  ,
     I         'CON'   , '   ' , 0     , 0     , 'X'     ,'(F12.3)' )
*
      CALL HISTRG
     I       ( 'SENS'  , 'sensible heat flux ' ,'W/m**2', 'ASFC',
     I          VMISS  , VMISS , 20.   , 100.  , ISTYPL  ,
     I         'CON'   , '   ' , 0     , 0     , 'X'     ,'(F12.3)' )
*
      CALL HISTRG
     I       ( 'QRADL' , 'longwave heating   ' ,'K/s    ', 'ALEV',
     I          VMISS  , VMISS , VMISS , VMISS , ISTYPL  ,
     I         'CON'   , '   ' , 0     , 0     , 'X'     ,'(1PE12.3)' )
*
      CALL HISTRG
     I       ( 'OLR'   , 'outgoing longwave  ' ,'W/m**2', 'ASFC',
     I          VMISS  , VMISS , VMISS , VMISS , ISTYPL  ,
     I         'CON'   , '   ' , 0     , 0     , 'X'     ,'(F12.3)' )
*
      CALL HISTRG
     I       ( 'TRPPM' , 'trpopause level' ,'mb    ', 'ASFC',
     I          VMISS  , VMISS , VMISS , VMISS , ISTYPL  ,
     I         'CON'   , '   ' , 0     , 0     , 'X'     ,'(1PE12.3)' )
*
      CALL HISTRG
     I       ( 'SLR'   , 'surface longwave   ' ,'W/m**2', 'ASFC',
     I          VMISS  , VMISS , VMISS , VMISS , ISTYPL  ,
     I         'CON'   , '   ' , 0     , 0     , 'X'     ,'(F12.3)' )
*
      CALL HISTRG
     I       ( 'QCND'  , 'condensation heating','K/s   ', 'ALEV',
     I          VMISS  , VMISS , VMISS , VMISS , ISTYPL  ,
     I         'CON'   , '   ' , 0     , 0     , 'X'     ,'(1PE12.3)' )
*
      CALL HISTRG
     I       ( 'SCND'  , 'condensation moistning','1/s   ', 'ALEV',
     I          VMISS  , VMISS , VMISS , VMISS , ISTYPL  ,
     I         'CON'   , '   ' , 0     , 0     , 'X'     ,'(1PE12.3)' )
*
      CALL HISTRG
     I       ( 'QCNDC' , 'cumulus heating'     ,'K/s   ', 'ALEV',
     I          VMISS  , VMISS , VMISS , VMISS , ISTYPL  ,
     I         'CON'   , '   ' , 0     , 0     , 'X'     ,'(1PE12.3)' )
*
      CALL HISTRG
     I       ( 'SCNDC' , 'cumulus moistning'   ,'1/s   ', 'ALEV',
     I          VMISS  , VMISS , VMISS , VMISS , ISTYPL  ,
     I         'CON'   , '   ' , 0     , 0     , 'X'     ,'(1PE12.3)' )
*
      CALL HISTRG
     I       ( 'SRMNQ' , 'remove negative moist','1/s   ', 'ALEV',
     I          VMISS  , VMISS , VMISS , VMISS , ISTYPL  ,
     I         'CON'   , '   ' , 0     , 0     , 'X'     ,'(1PE12.3)' )
*
      CALL HISTRG                              !" 保坂 93/08/29 追加
     I       ( 'LRAD' , 'LONG RADIATION' ,'W/m**2  '     , 'AMLEV',
     I          VMISS  , VMISS , VMISS , VMISS , ISTYPL  ,
     I         'CON'   , '   ' , 0     , 0     , 'X'     ,'(1PE12.3)' )
*
      CALL HISTRG                              !" 保坂 93/08/29 追加
     I       ( 'RUFLUX' , 'UPWARD RAD FLUX' ,'W/m**2  '     , 'AMLEV',
     I          VMISS  , VMISS , VMISS , VMISS , ISTYPL  ,
     I         'CON'   , '   ' , 0     , 0     , 'X'     ,'(1PE12.3)' )
*
      CALL HISTRG                              !" 保坂 93/08/29 追加
     I       ( 'RDFLUX' , 'DOWNWARD RAD FLUX' ,'W/m**2  '     , 'AMLEV',
     I          VMISS  , VMISS , VMISS , VMISS , ISTYPL  ,
     I         'CON'   , '   ' , 0     , 0     , 'X'     ,'(1PE12.3)' )
*
      CALL HISTRG                              !" 保坂 93/08/29 追加
     I       ( 'ALRAD' , 'LONG RADIATION' ,'W/m**2  '     , 'AMLEV',
     I          VMISS  , VMISS , VMISS , VMISS , ISTYPL  ,
     I         'CON'   , '   ' , 0     , 0     , 'X'     ,'(1PE12.3)' )
*
      CALL HISTRG                              !" 保坂 93/08/29 追加
     I       ( 'AUFLUX' , 'UPWARD RAD FLUX' ,'W/m**2  '     , 'AMLEV',
     I          VMISS  , VMISS , VMISS , VMISS , ISTYPL  ,
     I         'CON'   , '   ' , 0     , 0     , 'X'     ,'(1PE12.3)' )
*
      CALL HISTRG                              !" 保坂 93/08/29 追加
     I       ( 'ADFLUX' , 'DOWNWARD RAD FLUX' ,'W/m**2  '     , 'AMLEV',
     I          VMISS  , VMISS , VMISS , VMISS , ISTYPL  ,
     I         'CON'   , '   ' , 0     , 0     , 'X'     ,'(1PE12.3)' )
*
      RETURN
      END
