!================================================================
! 一次元放射対流平衡大気計算; メインプログラム
!
! 97/02/02 小高正嗣
! 97/02/18 小高正嗣; 地表温度を与える版
! 97/02/19 小高正嗣; 湿潤大気へ
! 97/03/12 小高正嗣; バグフィックス
! 97/03/16 小高正嗣; Ts を変数に 
! 97/05/04 小高正嗣; スキームテスト(p,tau,T,F は格子点上)
! 97/06/06 小高正嗣; 乾燥大気
!================================================================
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 :: xv	! 
	REAL,DIMENSION(:),ALLOCATABLE :: tau	! 光学的厚さ
	REAL			      :: dpln	! 圧力刻み(log)
	REAL			      :: Ftop	! 大気上端上向き放射束密度

!-内部変数

	REAL			      :: slp	! 全地表気圧
	REAL			      :: Ts	! 地表気温
	INTEGER			      :: nnl	! 鉛直層数
	INTEGER			      :: i

	dpln = LOG(10.)/100			! 圧力刻みの設定

	WRITE(*,*) 'INPUT SURFACE TEMP.(R) ?;'
	READ(*,*) Ts				!-地表気温読み込み

	slp = Ps 				!-地表全圧設定
	nnl = INT((LOG(slp) - LOG(Pt))/dpln) + 1!-鉛直層数設定

	WRITE(*,*) 'dpln=',dpln
	WRITE(*,*) 'nnl=',nnl

!-配列の割り付け: 第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) )

	CALL set_dprs( dpln, slp, p )		!-乾燥大気鉛直差分
	CALL dadiabat( p, Ts, ptemp )		!-乾燥対流圏温位計算
	CALL set_dtau( p, tau )			!-光学的厚さの計算(乾燥)
	CALL get_urf( tau, ptemp, Fu )		!-各層の上向き放射
	CALL get_drf( tau, ptemp, Fd )		!-各層の下向き放射
	CALL get_fnet( Fu, Fd, Fnet )		!-各層の正味上向き放射束密度

	DO i = 1,nnl

	   LTRP = i				!-圏界面高度の仮定
	   temp = ptemp
	   CALL rdeqiv( LTRP, tau, bf, temp, Ftop)!-成層圏の温度計算; 放射平衡
	   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 )	!-各層の正味上向き放射束密度

	   IF ( Fcnv(LTRP) < 0. ) THEN		!-対流熱流束=0を圏界面とする
		LTRP = i-1
		CALL rdeqiv( LTRP, tau, bf, temp, Ftop)
		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 )
		EXIT
	   END IF

	END DO

	WRITE(*,*) 'LTRP    =',LTRP
	WRITE(*,*) 'Ptp   =',p(LTRP)
	WRITE(*,*) 'tautp   =',tau(LTRP)
	WRITE(*,*) 'Ftop  =',Ftop
	WRITE(*,*) 'Btop  =',bf(nnl)
	WRITE(*,*) 'FUtop =',Fu(nnl)
	WRITE(*,*) 'FDtop =',Fd(nnl)
	WRITE(*,*) 'FCNVtp  =',Fcnv(LTRP)

	CALL dclfig( p,tau,bf,Fu,Fd,Fsol,Fcnv,temp)!-DCLでお絵書き

	WRITE(*,*) 'END PROGRAM'

	CONTAINS
!----------------------------------------------------------------------
SUBROUTINE dclfig(p,tau,bf,Fu,Fd,Fsol,Fcnv,temp)

	REAL,DIMENSION(:),INTENT(in) :: Fu	! 上向き放射束密度
	REAL,DIMENSION(:),INTENT(in) :: Fd	! 下向き放射束密度
	REAL,DIMENSION(:),INTENT(in) :: bf	! 黒体放射束密度
	REAL,DIMENSION(:),INTENT(in) :: Fsol	! 太陽放射束密度
	REAL,DIMENSION(:),INTENT(in) :: Fcnv	! 対流熱流束密度
	REAL,DIMENSION(:),INTENT(in) :: temp	! 温度
	REAL,DIMENSION(:),INTENT(in) :: p	! 圧力
	REAL,DIMENSION(:),INTENT(in) :: tau	! 光学的深さ

	INTEGER :: zdim
	INTEGER :: iws 				! 電脳パラメータ
	REAL :: taub 				! 光学的深さ
	REAL :: fnetmax
	REAL :: tempmax

	zdim = SIZE(p)-2
	fnetmax = MAXVAL(Fnet)
	tempmax = MAXVAL(temp)
	
	taub = tau(1)
	
	CALL swiset( 'IWIDTH', 700 )
	CALL swiset( 'IHEIGHT', 400 )
        CALL sglset( 'LCORNER', .FALSE. )
	WRITE(*,*) 'WORKSTATION ID (I) ?;'
	CALL sgpwsn
	READ(*,*) iws
	CALL gropn( -iws )
	CALL sldiv( 'Y',1, 2 )

	CALL grfrm
	CALL grswnd( 0., 600., taub, tau(zdim) )
	CALL grsvpt( 0.2, 0.8, 0.2, 0.8 )
	CALL grstrn( 2 )
	CALL grstrf
	CALL usdaxs
	CALL uxsttl( 'B', 'netflux (W/m^2)', 0. )
	CALL uysttl( 'L', 'optical depth', 0. )
	CALL uuslni(2)
	CALL uulin( zdim, Fu - Fd, tau )
	CALL uuslni(1)
	CALL uulin( zdim, Fu, tau )
	CALL uuslnt(2)
	CALL uulin( zdim, Fd, tau )
	CALL uuslnt(3)
	CALL uulin( zdim, Fcnv, tau )
	CALL uuslnt(1)
	CALL uuslni(2)

	CALL grfrm
	CALL grswnd( 200., tempmax+50., slp, p(zdim) )
	CALL grsvpt( 0.2, 0.8, 0.2, 0.8 )
	CALL grstrn( 2 )
	CALL grstrf
	CALL usdaxs
	CALL uxsttl( 'B', 'temperature (K)', 0. )
	CALL uysttl( 'L', 'pressure (Pa)', 0. )
	CALL uulin( zdim, temp, p )

	CALL grcls

END SUBROUTINE dclfig
!----------------------------------------------------------------------
END PROGRAM main
