!================================================================
! 一次元放射対流平衡大気計算; メインプログラム
!
! 97/02/02 小高正嗣
! 97/02/18 小高正嗣; 地表温度を与える版
! 97/02/19 小高正嗣; 湿潤大気へ
! 97/03/12 小高正嗣; バグフィックス
! 97/03/16 小高正嗣; Ts を変数に 
! 97/05/04 小高正嗣; スキームテスト(p,tau,T,F は格子点上)
! 97/05/31 小高正嗣; スキームテスト,各Tsに対するIRの計算
! 97/06/22 小高正嗣; 
! 97/06/23 小高正嗣; Ps を可変に
!================================================================
PROGRAM main
USE para_module
USE cod_module
USE rad_module
USE adjust_module
IMPLICIT NONE

!-共通変数

	INTEGER			      :: LTRP	! 圏界面高度番号
	REAL,DIMENSION(:),ALLOCATABLE :: Fu	! 上向き放射束密度
	REAL,DIMENSION(:),ALLOCATABLE :: Fd	! 下向き放射束密度
	REAL,DIMENSION(:),ALLOCATABLE :: bf	! 黒体放射束密度
	REAL,DIMENSION(:),ALLOCATABLE :: Fcnv	! 対流熱流束密度
	REAL,DIMENSION(:),ALLOCATABLE :: Fnet	! 正味上向き放射束密度
	REAL,DIMENSION(:),ALLOCATABLE :: Fsol	! 太陽放射束密度
	REAL,DIMENSION(:),ALLOCATABLE :: temp	! 温度
	REAL,DIMENSION(:),ALLOCATABLE :: ptemp	! 温位
	REAL,DIMENSION(:),ALLOCATABLE :: p	! 全圧力
	REAL,DIMENSION(:),ALLOCATABLE :: tau	! 光学的厚さ
	REAL,DIMENSION(:),ALLOCATABLE :: xv	! 水蒸気モル分率
	REAL,DIMENSION(:),ALLOCATABLE :: pxv	! 水蒸気モル分率
	REAL			      :: dpln	! 圧力刻み(log)
	REAL			      :: Ftop	! 大気上端上向き放射束密度

!-内部変数

	REAL			      :: slp	! 全地表気圧
	REAL			      :: PBR
	INTEGER			      :: nnl	! 鉛直層数
	INTEGER			      :: i,j,k
	INTEGER,PARAMETER	      :: dim = 31
	REAL,DIMENSION(5)	      :: Ps	! 非凝結成分地表気圧(Pa) 
	REAL,DIMENSION(dim)	      :: IR	! 大気上端正味放射量
	REAL,DIMENSION(dim)	      :: Ts	! 地表気温(配列)
	REAL,DIMENSION(dim)	      :: TBB	! 等価黒対輻射
	REAL 			      :: Tmin
	REAL			      :: Tmax
	REAL			      :: IRmin
	REAL			      :: IRmax

	dpln = LOG(10.)/100			! 圧力刻みの設定
	Tmin = 250.
	Tmax = 550.
	IRmin= 200.
	IRmax= 400.
	Ps(1)=0. ;Ps(2)=1.e+4 ;Ps(3)=1.e+5 ;Ps(4)=1.e+6 ;Ps(5)=1.e+7

	CALL dclopn				!-DCLの設定

	DO i = 1,5

	   WRITE(*,*) 'Ps    =',Ps(i)

	   DO j = 25,55,1			!-地表気温は 250-550K,10K刻み

	      Ts(j-24) = j*10.			!-地表気温設定
	      slp = Ps(i) + pv00*EXP(-ql/(R*Ts(j-24)))
						!-地表全圧設定
	      nnl = INT((LOG(slp) - LOG(Pt))/dpln) + 1
						!-鉛直層数設定

!-配列の割り付け: 第0層と第(nnl+1)層は境界条件のため
!		  フラックスは第1層から第(nnl)層まで

              ALLOCATE( Fu(nnl) )
	      ALLOCATE( Fd(nnl) )
	      ALLOCATE( bf(0:nnl+1) )
	      ALLOCATE( Fcnv(nnl) )
	      ALLOCATE( Fnet(0:nnl+1) )
	      ALLOCATE( Fsol(nnl) )
	      ALLOCATE( temp(0:nnl+1) )
	      ALLOCATE( ptemp(0:nnl+1) )
	      ALLOCATE( p(0:nnl+1) )
	      ALLOCATE( tau(0:nnl+1) )
	      ALLOCATE( xv(0:nnl+1) )
	      ALLOCATE( pxv(0:nnl+1) )

	      CALL set_mprs( dpln, slp, p )	!-湿潤大気鉛直差分
	      CALL madiabat( p, Ts(j-24), ptemp, pxv )
						!-湿潤対流圏温位計算
	      xv = pxv

	      CALL set_mtau( p, xv, tau )	!-光学的厚さの計算(湿潤)
	      CALL get_urf( tau, ptemp, Fu )	!-各層の上向き放射
	      CALL get_drf( tau, ptemp, Fd )	!-各層の下向き放射
	      CALL get_fnet( Fu, Fd, Fnet )	!-各層の正味上向き放射束密度

	      LTRP = nnl 

	      CALL fnetmaxl( Fnet, LTRP )	!-対流圏界面の決定;
						! 正味放射最大高度
	      CALL get_bbf( ptemp, bf )		!-

	      DO k = LTRP,1,-1	

	         PBR = ( Fu(LTRP) - bf(LTRP) )*(1.5*tau(LTRP) + 1.) &
	              - bf(LTRP)

	         IF ( PBR > 0. ) THEN
	            LTRP = LTRP - 1
 	         ELSE
	            EXIT
	         END IF

	      END DO

	      temp(0:LTRP) = ptemp(0:LTRP)	!-温度の設定; 断熱温度

	      CALL set_xv( LTRP, pxv, xv )	!-成層圏水蒸気モル分率の計算
	      CALL set_mtau( p, xv, tau )	!-光学的厚さの計算(湿潤)
	      CALL rdeqiv( LTRP, tau, bf, temp, Ftop, Fu, Fd )	
						!-成層圏の温度計算; 放射平衡
	      CALL get_urf( tau, temp, Fu )	!-各層の上向き放射
	      CALL get_drf( tau, temp, Fd )	!-各層の下向き放射
	      CALL get_bbf( temp, bf )		!-

 	      Fsol = Fu(nnl) - Fd(nnl)

	      CALL get_cvf( Fu, Fd, Fsol, Fcnv )	!-各層の対流熱流束
	      CALL get_fnet( Fu, Fd, Fnet )	!-各層の正味上向き放射束密度

	      IR(j-24) = Fnet(nnl)

 	      DEALLOCATE(Fu,Fd,bf,Fsol,Fcnv,Fnet,temp,ptemp,p,tau,xv,pxv)

	   END DO
	   CALL uulin( dim, Ts, IR )
	END DO

	CALL get_bbf( Ts, TBB )
	CALL uulin( 5, Ts, TBB )
	CALL grcls

	WRITE(*,*) 'END PROGRAM'

	CONTAINS
!----------------------------------------------------------------------
SUBROUTINE dclopn

	INTEGER :: iws 					! 電脳パラメータ

	CALL swiset( 'IWIDTH', 700 )
	CALL swiset( 'IHEIGHT', 400 )
        CALL sglset( 'LCORNER', .FALSE. )
	WRITE(*,*) 'WORKSTATION ID (I) ?;'
	CALL sgpwsn
	READ(*,*) iws

	CALL gropn( iws )
	CALL grfrm
	CALL grswnd( Tmin, Tmax, IRmin, IRmax )
	CALL grsvpt( 0.15, 0.85, 0.2, 0.8 )
	CALL grstrn( 1 )
	CALL grstrf
	CALL usdaxs
	CALL uxsttl( 'B', 'Surface Temp.(K)', 0. )
	CALL uysttl( 'L', 'IR (W/m2)', 0. )

END SUBROUTINE dclopn
!----------------------------------------------------------------------
END PROGRAM main
