!=================================================================
!  鉛直 1 次元火星大気放射対流平衡モデル: メインプログラム
!
!  1999/10/17 小高正嗣
!=================================================================
      program mars1D

#include "grid_size_M.f"
!-----------------------------------------------------------------      
      implicit real*8 ( a-h,o-z )
!-----------------------------------------------------------------      

      dimension
     \   PRES0(-2:NZB), TEMP0(-2:NZB), TPOT0(-2:NZB),
     \   PPAI0(-2:NZB), DENS0(-2:NZB), 
     \   RVAP0(-2:NZB), QVAP0(-2:NZB),
     \   TPOT1(-2:NZB)

      dimension
     \   DTPRAD(-2:NZB),
     \   DTPDI0(-2:NZB)

      dimension
     \   BBRAD ( -2:NZB ),
     \   FXRADU( -2:NZB ),
     \   FXRADD( -2:NZB ),
     \   FXRADN( -2:NZB )

      dimension 
     \   TKENRG(-2:NZB),
     \   CDTURB(-2:NZB), 
     \   CD1   (-2:NZB), 
     \   CD2   (-2:NZB),
     \   STABZ (-2:NZB),
     \   SSTURB(-2:NZB),
     \   DTQDI0(-2:NZB)

      dimension
     \   FHEAT ( 2 )

      dimension
     \   TSFC( -2:NXB ), QVSFC( -2:NXB ), FCOLI( -2:NXB )

      dimension 
     \   ZGRID(-2:NZB), FDZ(-2:NZB), FDZM(-2:NZB), OPL(-2:NZB)

      dimension EQWDTH(0:NZ,0:NZ)
      dimension work(-2:NZB)

!-----------------------------------------------------------------      
! 物理定数の設定
!-----------------------------------------------------------------      

      call SETCST
     O   ( GRAV, PSD, AMDRY, RDRY, CPDRY, RMVAP, ALATNT )

      call SETCS2
     O   ( CDRAG, VSFC0, TSFC, QVSFC,
     O     FCOLI,
     O     PSFC,
     C     RMVAP,
     C     FXSOL,
     D     NX, NXB                      )

!-----------------------------------------------------------------      
! 鉛直格子の設定
!-----------------------------------------------------------------      

      call SETGRD
     O   ( ZGRID, FDZ, FDZM, DX,
     D     NZ, NZB                )

!-----------------------------------------------------------------      
! 初期温度分布の設定
!-----------------------------------------------------------------      

      do IZ = 0, NZ
         TEMP0( IZ ) = TEMPZ( ZGRID(IZ) )
      end do

      call BOUND ( TEMP0, NZ, NZB )

!-----------------------------------------------------------------      
! 初期拡散係数の設定
!-----------------------------------------------------------------      

      do IZ = 0, NZ
         CDTURB( IZ ) = 1.0D0
      end do

      call BOUND ( CDTURB, NZ, NZB )

!-----------------------------------------------------------------      
! ループ計算開始
!-----------------------------------------------------------------      

      DTI = 10.0D0

c      do 

         call VPROF             ! 大気構造の計算
     I      ( ZGRID, PSFC,
     O        PRES0, DENS0, TEMP0, TPOT0, PPAI0,
     O        QVAP0, RVAP0,
     C        GRAV , RDRY , CPDRY, PSD,  RMVAP,
     G        FDZ  , FDZM , DX   ,
     D        NZ   , NZB                         )


         call SETOPL            ! 光路長の計算
     I      ( PRES0, PSFC  ,
     O        OPL  , 
     C        GRAV , 
     D        NZ   , NZB     )

         call CLEQWD            ! 等価幅の計算
     I      ( TEMP0 , PRES0, OPL , 
     I        ZGRID , FDZ  , FDZM, DX, 
     O        EQWDTH,                       
     D        NZ    , NZB      )

         call CDCOEF            ! 半整数格子上の拡散係数の計算(含リミッター)
     I      ( CDTURB,
     O        CD1   , CD2   ,
     B        ZGRID ,
     D        NZ    , NZB     )

         call CLRAD             ! 放射加熱率の計算
     I      ( TPOT0 , EQWDTH, PSFC,
     O        DTPRAD, FXRADU, FXRADD, FXRADN, BBRAD,
     B        TEMP0 , PRES0 , DENS0, FDZ  ,
     C        TSFC  , CPDRY , RDRY , DX   , 
     D        NX    , NXB   , NZ   , NZB     )

         call ADDIF0            ! 拡散加熱率の計算
     I      ( TPOT0 ,  CD1 ,
     O        DTPDI0,
     B        DENS0,
     G        DX, FDZ, FDZM,
     D        NZ, NZB        )

         call FXSURF            ! 地表熱フラックスの計算
     I      ( TPOT , CD2  ,
     O        FHEAT, 
     B        DENS0, TPOT0, PPAI0,
     C        CPDRY,
     S        TSFC , 
     D        NX   , NXB  , NZ   , NZB )

         do IZ = 0, NZ-1
            TPOT1(IZ) = TPOT0(IZ)
     \           + DTI  * ( DTPDI0(IZ) + DTPRAD(IZ)  )
         end do
       
         TPOT1(0) = TPOT1(0)
     \     + DTI * FHEAT(1) / DX / FDZ(0) / DENS0(0) / CPDRY
         TPOT1(NZ-1) = TPOT1(NZ-1)
     \     + DTI * FHEAT(2) / DX / FDZ(NZ-1) / DENS0(NZ-1) / CPDRY

         do IZ = 0, NZ-1
            TPOT0( IZ ) = TPOT1( IZ )
            TEMP0( IZ ) = TPOT0( IZ ) * PPAI0( IZ )
         end do

         call BOUND( TPOT0, NZ, NZB )
         call BOUND( TEMP0, NZ, NZB )


!-----------------------------------------------------------------      
! 地表面温度の計算
!-----------------------------------------------------------------      

         call CLTSFC
     I   ( FHEAT , FXRADD, FXSOL,
     O     TSFC  ,
     C     SIGMA ,
     D     NX    , NXB   , 
     D     NZ    , NZB            )

!-----------------------------------------------------------------      
! 乱流拡散係数の計算
!-----------------------------------------------------------------      

         call ADDIF0            ! 拡散加熱率の計算
     I      ( TKENRG,  CD2 ,
     O        DTQDI0,
     B        DENS0,
     G        DX, FDZ, FDZM,
     D        NZ, NZB )

         call CLTURB            ! 乱流生成と散逸の計算
     O      ( SSTURB,
     I        TKENRG, CDTURB,
     B        TPOT0 , 
     G        ZGRID , FDZ   , FDZM  , DX    ,
     C        GRAV  , 
     W        STABZ , 
     D        NX    , NZ    , NXB   , NZB    )


         do IZ = 0, NZ-1
            TKENRG( IZ ) = TKENRG( IZ ) 
     \           + DTI * ( DTQDI0(IZ) + SSTURB(IZ) )
         end do

         do IZ = 0, NZ-1

            MIXLNG = MIN( ZGRID(IZ), FDZ(IZ) )

            CDTURB( IZ ) = 
     \           MAX( 0.0D0, MIN( 800.0D0, 
     \                0.2D0 * sqrt( TKENRG(IZ) ) *  MIXLNG ) )

            TKENRG( IZ ) = ( CDTURB( IZ ) / 0.2D0 / MIXLNG )**2.0D0

         end do

         call BOUND( TKENRG, NZ, NZB )
         call BOUND( CDTURB, NZ, NZB )

c      end do

         do iz = -2, NZB
            DTPRAD(iz) = DTPRAD(iz) * 86400.0D0
         end do

         do iz = 0,NZ
            WORK(iz) = EQWDTH(iz,50)
         end do

         call dclfig
     I     (  ZGRID , WORK,
     D        NZ    , NZB            )

         call dclfig
     I     (  ZGRID , FXRADU,
     D        NZ    , NZB            )

         do iz = 198,202
            write(*,*) IZ, FDZ(IZ)
         end do

         do iz = 218,222
            write(*,*) IZ, FDZ(IZ)
         end do

c      write(*,*) 'IZ   ','ZGRID     ','FXRADU      ','FXRADD      ',
c     \        'FXRADN      ','Qrad     ', 'DENS0'
c      do iz = NZ-20, NZB
c      write(*,*) iz, ZGRID(iz), FXRADU(iz), FXRADD(iz), 
c     \        FXRADN(iz),DTPRAD(iz),DENS0(iz)
c      end do         

c         call dclfig
c     I   ( ZGRID , TEMP0 , PRES0 , OPL   ,
c     I     FXRADU, FXRADD, FXRADN, BBRAD ,
c     I     DTPRAD,
c     D     NZ    , NZB                     )


      end

