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

!-共通変数

	INTEGER              :: TP	! 圏界面高度番号
	REAL,DIMENSION(zdim) :: uf	! 上向き放射束密度
	REAL,DIMENSION(zdim) :: df	! 下向き放射束密度
	REAL,DIMENSION(zdim) :: bf	! 黒体放射束密度
	REAL,DIMENSION(zdim) :: sf	! 太陽放射束密度
	REAL,DIMENSION(zdim) :: cf	! 対流熱流束密度
	REAL,DIMENSION(zdim) :: temp	! 温度
	REAL,DIMENSION(zdim) :: ptemp	! 温位
	REAL,DIMENSION(zdim) :: p	! 全圧力
	REAL,DIMENSION(zdim) :: tau	! 光学的厚さ
	REAL,DIMENSION(zdim) :: xv	! 水蒸気モル分率
	REAL,DIMENSION(zdim) :: pxv	! 水蒸気モル分率
	REAL                 :: uf0	! 大気上端上向き放射束密度
	REAL                 :: slp	! 全地表気圧
	REAL                 :: Ts	! 地表気温
	INTEGER		     :: i

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

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

	DO i = 2,zdim

	  TP = i
	  temp = ptemp
	  CALL get_bbf(temp,bf)			!-黒体放射計算

	  IF ( TP == zdim ) THEN		!-全層が対流圏の場合
	    WRITE(*,*) 'ALL LAYER IS TROPOSPHERE'

	    TP = TP - 1
	    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_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(zdim)
	    WRITE(*,*) 'UFtop=',uf(zdim)
	    WRITE(*,*) 'DFtop=',df(zdim)
	    WRITE(*,*) 'CFtop=',cf(zdim)

	    EXIT	  
	  END IF

	  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)		!-圏界面での対流熱流束

	  IF ( cf(TP) <= 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(zdim)
	    WRITE(*,*) 'UFtop=',uf(zdim)
	    WRITE(*,*) 'DFtop=',df(zdim)
	    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(zdim),INTENT(in) :: uf	! 上向き放射束密度
	REAL,DIMENSION(zdim),INTENT(in) :: df	! 下向き放射束密度
	REAL,DIMENSION(zdim),INTENT(in) :: bf	! 黒体放射束密度
	REAL,DIMENSION(zdim),INTENT(in) :: sf	! 太陽放射束密度
	REAL,DIMENSION(zdim),INTENT(in) :: cf	! 対流熱流束密度
	REAL,DIMENSION(zdim),INTENT(in) :: temp	! 温度
	REAL,DIMENSION(zdim),INTENT(in) :: p	! 圧力
	REAL,DIMENSION(zdim),INTENT(in) :: tau	! 光学的深さ
	REAL,DIMENSION(zdim),INTENT(in) :: xv	! 水蒸気モル分率

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

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

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

	CALL grfrm
!	CALL grswnd( 0., temp(1)+100., taub, 0. )
	CALL grswnd( 0., temp(1)+100., taub, tau(zdim-1) )
	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-1, temp, tau )

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

	CALL grfrm
	CALL grswnd( 1.e-8, 1., taub, 0. )
	CALL grsvpt( 0.2, 0.8, 0.2, 0.8 )
	CALL grstrn( 3 )
	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-1) )
	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-1, xv, p )

	CALL grcls

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