!================================================================
! 一次元放射対流平衡大気計算; メインプログラム
!
! 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			      :: TP	! 圏界面高度番号
	REAL,DIMENSION(:),ALLOCATABLE :: uf	! 上向き放射束密度
	REAL,DIMENSION(:),ALLOCATABLE :: df	! 下向き放射束密度
	REAL,DIMENSION(:),ALLOCATABLE :: bf	! 黒体放射束密度
	REAL,DIMENSION(:),ALLOCATABLE :: sf	! 太陽放射束密度
	REAL,DIMENSION(:),ALLOCATABLE :: cf	! 対流熱流束密度
	REAL,DIMENSION(:),ALLOCATABLE :: Fnet	! 正味上向き放射束密度
	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			      :: uf0	! 大気上端上向き放射束密度
	REAL			      :: slp	! 全地表気圧
	REAL			      :: Ts	! 地表気温
	INTEGER			      :: nnl	! 鉛直層数
	INTEGER			      :: i

	dpln = LOG(10.)/100			! 圧力刻みの設定
	WRITE(*,*) 'dpln=',dpln
	WRITE(*,*) 'INPUT SURFACE TEMP.(R) ?;'
	READ(*,*) Ts				!-地表気温読み込み

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

	WRITE(*,*) 'nnl=',nnl
!-配列の割り付け

	ALLOCATE(uf(nnl))
	ALLOCATE(df(nnl))
	ALLOCATE(bf(nnl))
	ALLOCATE(sf(nnl))
	ALLOCATE(cf(nnl))
	ALLOCATE(Fnet(nnl))
	ALLOCATE(temp(nnl))
	ALLOCATE(ptemp(nnl))
	ALLOCATE(p(nnl))
	ALLOCATE(tau(nnl))
	ALLOCATE(xv(nnl))
	ALLOCATE(pxv(nnl))

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

	DO i = nnl-1,1,-1

	  TP = i
	  write(*,*) TP
	  temp = ptemp
	  CALL get_bbf(temp,bf)			!-黒体放射計算
	  CALL set_xv(TP,pxv,xv)		!-水蒸気モル分率の計算
	  CALL set_mtau(p,xv,tau)		!-光学的厚さの計算
!	  CALL set_dtau(p,tau)			!-光学的厚さの計算
	  CALL rdeqiv(TP,tau,bf,temp,uf0)	!-成層圏の温度計算; 放射平衡
	  sf = uf0
	  CALL get_urf1(TP,tau,temp,uf)		!-圏界面での上向き放射
	  CALL get_drf1(TP,tau,temp,df)		!-圏界面での下向き放射
	  CALL get_cvf(uf,df,sf,cf)		!-圏界面での対流熱流束
	  CALL get_fnet(uf,df,Fnet)		!-全層での正味上向き放射束密度

!	  IF ( cf(TP) <= 0.)  THEN		!-対流圏界面の決定
	  IF ( Fnet(TP)-Fnet(TP+1) <= 0.)  THEN	!-対流圏界面の決定
	    TP = TP + 1				!-一つ下が圏界面
	    CALL get_bbf(temp,bf)
	    CALL set_xv(TP,pxv,xv)
	    CALL set_mtau(p,xv,tau)
	    CALL rdeqiv(TP,tau,bf,temp,uf0)
	    sf = uf0
	    CALL get_urf(tau,temp,uf)		!-全層での上向き放射
	    CALL get_drf(tau,temp,df)		!-全層での下向き放射
	    CALL get_cvf(uf,df,sf,cf)		!-全層での対流熱流束

	    WRITE(*,*) 'TP   =',TP
	    WRITE(*,*) 'Ptp  =',p(TP)
	    WRITE(*,*) 'Ftop =',uf0
	    WRITE(*,*) 'Btop =',bf(nnl)
	    WRITE(*,*) 'UFtop=',uf(nnl)
	    WRITE(*,*) 'DFtop=',df(nnl)
!	    WRITE(*,*) 'CFtp=',cf(TP)

	    EXIT
	  END IF
	END DO

	CALL dclfig(p,tau,bf,uf,df,sf,cf,temp,xv)	!-DCLでお絵書き

	WRITE(*,*) 'END PROGRAM'

	CONTAINS
!----------------------------------------------------------------------
SUBROUTINE dclfig(p,tau,bf,uf,df,sf,cf,temp,xv)

	REAL,DIMENSION(:),INTENT(in) :: uf	! 上向き放射束密度
	REAL,DIMENSION(:),INTENT(in) :: df	! 下向き放射束密度
	REAL,DIMENSION(:),INTENT(in) :: bf	! 黒体放射束密度
	REAL,DIMENSION(:),INTENT(in) :: sf	! 太陽放射束密度
	REAL,DIMENSION(:),INTENT(in) :: cf	! 対流熱流束密度
	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)
	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., uf(1)+100., taub, 0. )
!	CALL grswnd( 0., uf(1)+100., taub, tau(nnl) )
	CALL grswnd( 0., uf(zdim)-df(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 uulin( zdim, uf - df, tau )
!	CALL uuslnt(2)
!	CALL uulin( zdim, df, tau )
!	CALL uuslnt(3)
!	CALL uulin( zdim, sf, tau )
!	CALL uuslnt(4)
!	CALL uulin( zdim, cf, tau )
!	CALL uuslnt(1)
!	CALL uuslni(2)
!	CALL uulin( zdimnnl, bf, tau )

	CALL grfrm
!	CALL grswnd( 0., uf(1)+100., slp, p(zdim) )
	CALL grswnd( 0., uf(zdim)-df(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 uuslnt(1)
	CALL uulin( zdim, uf - df, p )
!	CALL uuslnt(2)
!	CALL uulin( zdim, df, p )
!	CALL uuslnt(3)
!	CALL uulin( zdim, sf, p )
!	CALL uuslnt(4)
!	CALL uulin( zdim, cf, p )
!	CALL uuslnt(1)
!	CALL uuslni(2)
!	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(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, 1., 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, 1., 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
