!================================================================
! 一次元放射対流平衡大気計算; メインプログラム
!
! 97/02/02 小高正嗣
! 97/02/18 小高正嗣; 地表温度を与える版
! 97/02/19 小高正嗣; 湿潤大気へ
! 97/03/12 小高正嗣; バグフィックス
! 97/03/16 小高正嗣; Ts を変数に 
! 97/05/04 小高正嗣; スキームテスト(p,tau,T,F は格子点上)
!================================================================
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			      :: Ts	! 地表気温
	INTEGER			      :: nnl	! 鉛直層数

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

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

	slp = Ps + pv00*EXP(-ql/(R*Ts))		!-地表全圧設定
	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) )
	ALLOCATE( pxv(0:nnl+1) )

	CALL set_mprs( dpln, slp, p )		!-湿潤大気鉛直差分
	CALL madiabat( p, Ts, 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 )	!-対流圏界面の決定; 正味放射最大高度
	temp = ptemp				!-温度の設定; 断熱温度
	CALL set_xv( LTRP, pxv, xv )		!-成層圏水蒸気モル分率の計算
	CALL set_mtau( p, xv, tau )		!-光学的厚さの計算(湿潤)
	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 )		!-各層の正味上向き放射束密度

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

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

	WRITE(*,*) 'END PROGRAM'

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

	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	! 光学的深さ
	REAL,DIMENSION(:),INTENT(in) :: xv	! 水蒸気モル分率
	INTEGER :: zdim
	INTEGER :: iws 				! 電脳パラメータ
	REAL :: taub 				! 光学的深さ

	zdim = SIZE(p)-2

	taub = tau(1)
	
	CALL swiset( 'IWIDTH', 700 )
	CALL swiset( 'IHEIGHT', 800 )
        CALL sglset( 'LCORNER', .FALSE. )
	WRITE(*,*) 'WORKSTATION ID (I) ?;'
	CALL sgpwsn
	READ(*,*) iws
	CALL gropn( iws )
	CALL sldiv( 'Y',2, 3 )

	CALL grfrm
!	CALL grswnd( 0., Fu(1)+100., taub, 0. )
!	CALL grswnd( 0., Fu(1)+100., taub, tau(nnl) )
	CALL grswnd( 0., Fu(zdim)-Fd(zdim)+100., taub, tau(zdim) )
	CALL grsvpt( 0.2, 0.8, 0.2, 0.8 )
!	CALL grstrn( 1 )
	CALL grstrn( 2 )
	CALL grstrf
	CALL usdaxs
	CALL uxsttl( 'B', 'flux (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, Fsol, tau )
	CALL uuslnt(4)
	CALL uulin( zdim, bf, tau )
	CALL uuslnt(1)


	CALL grfrm
!	CALL grswnd( 0., Fu(1)+100., slp, p(zdim) )
	CALL grswnd( 0., Fu(zdim)-Fd(zdim)+100., slp, p(zdim) )
	CALL grsvpt( 0.2, 0.8, 0.2, 0.8 )
	CALL grstrn( 2 )
	CALL grstrf
	CALL usdaxs
	CALL uxsttl( 'B', 'flux (W/m^2)', 0. )
	CALL uysttl( 'L', 'pressure (Pa)', 0. )
	CALL uuslni(2)
	CALL uulin( zdim, Fu - Fd, p )
	CALL uuslni(1)
	CALL uulin( zdim, Fu, p )
	CALL uuslnt(2)
	CALL uulin( zdim, Fd, p )	
	CALL uuslnt(3)
	CALL uulin( zdim, Fsol, p )
	CALL uuslnt(4)
	CALL uulin( zdim, bf, p )

	CALL grfrm
!	CALL grswnd( 0., temp(1)+100., taub, 0. )
	CALL grswnd( 0., temp(1)+100., taub, tau(zdim) )
	CALL grsvpt( 0.2, 0.8, 0.2, 0.8 )
!	CALL grstrn( 1 )
	CALL grstrn( 2 )
	CALL grstrf
	CALL usdaxs
	CALL uxsttl( 'B', 'temperature (K)', 0. )
	CALL uysttl( 'L', 'optical depth', 0. )
	CALL uuslni(2)
	CALL uuslnt(1)
	CALL uulin( zdim, temp, tau )

	CALL grfrm
	CALL grswnd( 0., temp(1)+100., 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 uuslni(1)
	CALL uulin( zdim, temp, p )

	CALL grfrm
	CALL grswnd( 1.e-8, 10., taub, tau(zdim) )
	CALL grsvpt( 0.2, 0.8, 0.2, 0.8 )
	CALL grstrn( 4 )
	CALL grstrf
	CALL usdaxs
	CALL uxsttl( 'B', 'mol fraction', 0. )
	CALL uysttl( 'L', 'optical depth', 0. )
!	CALL uuslni(1)
	CALL uulin( zdim, xv, tau )

	CALL grfrm
	CALL grswnd( 1.e-8, 10., slp, p(zdim) )
	CALL grsvpt( 0.2, 0.8, 0.2, 0.8 )
	CALL grstrn( 4 )
	CALL grstrf
	CALL usdaxs
	CALL uxsttl( 'B', 'mol fraction', 0. )
	CALL uysttl( 'L', 'pressure (Pa)', 0. )
!	CALL uuslni(1)
	CALL uulin( zdim, xv, p )

	CALL grcls

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