
!BS1998 の放射プロファイルを生成
!2012/05/22

program BS1998rad

!暗黙の型宣言禁止

  implicit none


!変数定義

  real(8), allocatable :: Qsub(:)  ! 放射プロファイル
  real(8), allocatable :: z_Z(:)   ! 高度座標
  real(8) :: dz                    ! 格子間隔
  real(8) :: z_Upper = 6.70d4       ! 上部頂点高度
  real(8) :: z_Lower = 2.70d4       ! 下部頂点高度
  real(8) :: Zmax = 1.0d5       ! 上端高度
  real(8) :: Zmin = 0.0d0       ! 下端高度
  real(8) :: c_Upper = 2.70d-2     ! 第 1 係数
  real(8) :: c_Lower = 3.60d-3     ! 第 2 係数
  real(8) :: s_Upper = 7.50d3      ! 第 1 標準偏差
  real(8) :: s_Lower = 1.30d4      ! 第 2 標準偏差

  integer :: nz = 200 ! 格子点数
  integer :: iz, ins

! 配列割付
  
  allocate( Qsub(1:nz), z_Z(1:nz))


! 配列初期化

  Qsub = 0.0d0
  z_Z = 0.0d0
  z_Z(1) = Zmin

! 座標配列設定

  dz = ( Zmax - Zmin) / dble(nz)

  do iz = 2, nz
     z_Z(iz) = z_Z(iz-1) + dz
  enddo


! 放射プロファイル計算

  do iz = 1, nz
     Qsub(iz) = &
       & c_Lower * exp( - ( z_Z(iz) - z_Lower )**2.0d0 / (2.0d0 * s_Lower**2.0d0)) &
       & + c_Upper * exp( - ( z_Z(iz) - z_Upper )**2.0d0 / (2.0d0 * s_Upper**2.0d0))
  enddo

! DCL で描画

  WRITE(*,*) ' WORKSTATION ID (I)  ? ;'  
  CALL SGPWSN  
  READ (*,*) ins  
  CALL GROPN( ins )  
!    CALL GROPN( 2 )  

  CALL GRFRM  
  CALL USSTTL( 'Qsub', 'W.m-3', 'z-coordination', 'm' )  

  CALL USGRPH( nz, real(Qsub), real(z_Z) )  
  CALL GRCLS

deallocate( Qsub, z_Z )



end program BS1998rad
