!================================================================
! 一次元放射対流平衡大気計算; メインプログラム
!
! 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/21 小高正嗣; スキームテスト,各Tsに対するIRの計算
! 97/06/22 小高正嗣; スキームテスト,各Tsに対するTの計算
!================================================================
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	! 地表気温
	REAL			      :: PBR
	INTEGER,PARAMETER	      :: dim=11
	REAL,DIMENSION(dim)	      :: kaid	! 非凝結成分吸収係数(1/kg.m)
	INTEGER			      :: nnl	! 鉛直層数
	INTEGER			      :: i,j
	INTEGER,PARAMETER	      :: Tsmin=25
	INTEGER,PARAMETER	      :: Tsmax=35

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

!	CALL dclopn_tvsp(Tsmin, Tsmax)		!-DCL初期設定
	CALL dclopn_tvstau( Tsmin, Tsmax )
!	CALL dclopn_xvvsp(Tsmin, Tsmax)
!	CALL dclopn_xvvstau

	kaid(1)=1.e-6 ;kaid(2)=5.e-6 ;kaid(3)=1.e-5
	kaid(4)=5.e-5 ;kaid(5)=1.e-4 ;kaid(6)=5.e-4
	kaid(7)=1.e-3 ;kaid(8)=5.e-3 ;kaid(9)=1.e-2
	kaid(10)=5.e-2 ;kaid(11)=1.e-1

	DO i = 1,dim

	   OPEN(10,file='kaid.dat')
	   WRITE(10,*) kaid(i)
	   CLOSE(10)
	
	   Ts = 215.				!-地表気温設定
	   slp = Ps + pv00*EXP(-ql/(R*Ts))	!-地表全圧設定
	   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_dprs( dpln, slp, p )	!-湿潤大気鉛直差分
	   CALL dadiabat( slp, p, Ts, ptemp )	!-湿潤対流圏温位計算
!	   xv = pxv
	   CALL set_dtau( p, tau )		!-光学的厚さの計算(湿潤)
!	   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 j = 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 )	!-各層の正味上向き放射束密度

	   WRITE(*,*) 'kaid    =',kaid(i),',','Fnet  =',Fnet(nnl)

!	   CALL uulin( nnl+2, temp, p )		!-DCLでお絵書き
	   CALL uulin( nnl+2, temp, tau )	!-DCLでお絵書き
!	   CALL uulin( nnl+2, xv, p )		!-DCLでお絵書き
!	   CALL uulin( nnl+2, xv, tau )		!-DCLでお絵書き

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

	END DO

	CALL grcls

	WRITE(*,*) 'END PROGRAM'

	CONTAINS
!----------------------------------------------------------------------
SUBROUTINE dclopn_tvsp( min, max )

	INTEGER,INTENT(in)	:: min, max 
	INTEGER			:: iws		! 電脳パラメータ
	REAL :: Tmin,Tmax,pmax,pmin

	Tmin = 0.
	Tmax = max*10.+50.
	pmin = Pt
	pmax = Ps + pv00*EXP(-ql/(R*Tmax))

	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, pmax, pmin )
	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. )

END SUBROUTINE dclopn_tvsp
!----------------------------------------------------------------------
SUBROUTINE dclopn_tvstau( min, max )

	INTEGER,INTENT(in)	:: min, max 
	INTEGER			:: iws		! 電脳パラメータ
	REAL :: Tmin,Tmax,taumax,taumin

	Tmin = 0.
	Tmax = max*10.+50.
	taumin = 1.e-10
	taumax = 1.e+2

	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, taumax, taumin )
	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', 'optical depth', 0. )

END SUBROUTINE dclopn_tvstau
!----------------------------------------------------------------------
SUBROUTINE dclopn_xvvsp( min, max )

	INTEGER,INTENT(in)	:: min, max 
	INTEGER			:: iws		! 電脳パラメータ
	REAL :: xvmin,xvmax
	REAL :: Tmin,Tmax,pmax,pmin

	Tmin = 0.
	Tmax = max*10.+50.
	pmin = Pt
	pmax = Ps + pv00*EXP(-ql/(R*Tmax))
	xvmin = 1.e-8
	xvmax = 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 grfrm
	CALL grswnd( xvmin, xvmax, pmax, pmin )
	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. )


END SUBROUTINE dclopn_xvvsp
!----------------------------------------------------------------------
SUBROUTINE dclopn_xvvstau

	INTEGER			:: iws		! 電脳パラメータ
	REAL :: xvmin,xvmax,taumax,taumin

	xvmin = 1.e-8
	xvmax = 1.
	taumin = 1.e-10
	taumax = 1.e+2

	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( xvmin, xvmax, taumax, taumin )
	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. )

END SUBROUTINE dclopn_xvvstau
!----------------------------------------------------------------------
END PROGRAM main



