!=================================================================
!  一次元 CO2 大気放射対流平衡計算プログラム; 
!  等価幅の計算(Pollack et al., 1981:J. Atmos. Sci., 38, 3--29.)
!
!  1999/06/04  小高正嗣
!  2000/02/21  小高正嗣; 補正ルーチンを加える
!  2000/02/22  小高正嗣; 大気上端からの寄与を考慮
!  2000/05/01  小高正嗣; 補正ルーチンを変更
!=================================================================

      subroutine CLEQWD
     I   ( TEMP0 , PRES0, OPL , 
     I     ZGRID , FDZ  , FDZM, DX, 
     O     EQWDTH,
     D     NZ  , NZB      )

!-----------------------------------------------------------------      
      implicit real*8 ( a-h,o-z )
!-----------------------------------------------------------------      

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

      dimension 
     \   EQWDTH(0:NZ+1,0:NZ+1)

      parameter ( EPSILON = 1.0D-15 )
      parameter ( TREF    = 200.0D0 )
!      parameter ( PREF    = 700.0D0 )

!-----------------------------------------------------------------      
! 参照温度の計算;2000/05/01
!-----------------------------------------------------------------      

      T0MEAN = 0.0D0
      P0MEAN = 0.0D0

      ZLENG  = ( ZGRID( NZ-1 ) - FDZ( NZ-1 )* DX /2.0D0 ) -
     \     ( ZGRID( 0 ) - FDZ( 0 )* DX /2.0D0 ) 
            
      do iz = 0, NZ - 1
         T0MEAN = T0MEAN + 
     \        TEMP0( iz ) * FDZ( iz ) * DX / MAX(ZLENG,EPSILON)
         P0MEAN = P0MEAN + 
     \        PRES0( iz ) * FDZ( iz ) * DX / MAX(ZLENG,EPSILON)
      end do

!      TREF = T0MEAN
      PREF = P0MEAN

      write(*,*) T0MEAN, P0MEAN

!-----------------------------------------------------------------      

      T0MEAN = 200.0D0

      do iz = 0, NZ
         EQWDTH(iz,iz)=0.0D0
      end do
      
      do iz = 1, NZ-1

         do jz = 0, iz-1

            T0MEAN = 0.0D0
            P0MEAN = 0.0D0

            ZLENG  = ( ZGRID( iz ) - FDZ( iz )* DX /2.0D0 ) -
     \               ( ZGRID( jz ) - FDZ( jz )* DX /2.0D0 ) 
            
            do kz = jz, iz - 1
               T0MEAN = T0MEAN + 
     \              TEMP0( kz ) * FDZ( kz ) * DX / MAX(ZLENG,EPSILON)
               P0MEAN = P0MEAN + 
     \              PRES0( kz ) * FDZ( kz ) * DX / MAX(ZLENG,EPSILON)
            end do
           
!            T0MEAN = MAX( T0MEAN, TREF )
!            P0MEAN = MAX( P0MEAN, PREF )

            EQWDTH(iz,jz) = 54.4D0 * ( T0MEAN / TREF )**(0.879D0) *
     \           LOG( 1.0D0 + 0.153D0 * ( T0MEAN / TREF )**(-0.256D0) *
!     \                ( ABS ( OPL( jz ) - OPL( iz ) ) )**(0.566D0) *
     \      ( ABS ( OPL( jz ) - OPL( iz ) )/1.96*100.0D0 )**(0.566D0) *
     \                P0MEAN**(0.323D0) )

         end do
      end do

! モデル最上層との等価幅: ここでは一つ下と同じにする. 

      do iz = 0, NZ-1
         EQWDTH(NZ,iz) = EQWDTH(NZ-1,iz)
      end do

! 上三角成分の計算(成分は対称)

      do iz = 0, NZ
         do jz = iz+1, NZ
            EQWDTH(iz, jz) = EQWDTH(jz, iz)
         end do
      end do

!-----------------------------------------------------------------      
! 大気上端との間の等価幅の計算;2000/02/22
!-----------------------------------------------------------------      

      do jz = 0, NZ-1
         
         T0MEAN = 0.0D0
         P0MEAN = 0.0D0
         
         ZLENG  = ( ZGRID( NZ ) + FDZ( NZ )* DX /2.0D0 ) -
     \        ( ZGRID( jz ) - FDZ( jz )* DX /2.0D0 ) 
         
         do kz = jz, NZ-1
            T0MEAN = T0MEAN + 
     \           TEMP0( kz ) * FDZ( kz ) * DX / MAX(ZLENG,EPSILON)
            P0MEAN = P0MEAN + 
     \           PRES0( kz ) * FDZ( kz ) * DX / MAX(ZLENG,EPSILON)
         end do

!         T0MEAN = MAX( T0MEAN, TREF )
!         P0MEAN = MAX( P0MEAN, PREF )
            
         EQWDTH(NZ+1,jz) = 54.4D0 * ( T0MEAN / TREF )**(0.879D0) *
     \        LOG( 1.0D0 + 0.153D0 * ( T0MEAN / TREF )**(-0.256D0) *
!     \             ( ABS ( OPL( jz ) ) )**(0.566D0) *
     \             ( ABS ( OPL( jz ) ) / 1.96 * 1.0D+2 )**(0.566D0) *
     \        P0MEAN**(0.323D0) )
            
      end do

      EQWDTH(NZ+1,NZ) = 54.4D0 * ( T0MEAN / TREF )**(0.879D0) *
     \     LOG( 1.0D0 + 0.153D0 * ( T0MEAN / TREF )**(-0.256D0) *
!     \     ( ABS ( OPL( NZ ) ) )**(0.566D0) *
     \     ( ABS ( OPL( NZ ) ) / 1.96 )**(0.566D0) *
     \     P0MEAN**(0.323D0) )
            
      do iz = 0, NZ
         EQWDTH(iz, NZ+1) = EQWDTH(NZ+1, iz)
      end do

!-----------------------------------------------------------------      
! 上層での等価幅の修正;2000/02/21,2000/02/22
! ある高度より高いところで等価幅が減少してしまう問題点の補正
!-----------------------------------------------------------------      

      do iz = 0, NZ
         do jz = iz+1, NZ+1
            EQWDTH(iz,jz) = MAX( EQWDTH(iz,jz), EQWDTH(iz,jz-1) )
         end do
      end do

      do iz = 0, NZ
         do jz = iz+1, NZ+1
            EQWDTH(jz, iz) = EQWDTH(iz, jz)
         end do
      end do

      return

      end
