!================================================================
! 一次元放射対流平衡大気計算; 放射計算
!
! 97/02/23 小高正嗣
! 97/03/12 小高正嗣; 修正
! 97/03/17 小高正嗣
! 97/05/08 小高正嗣
!================================================================
MODULE rad_module
USE para_module
IMPLICIT NONE
PUBLIC :: get_bbf,get_urf1,get_drf1,get_urf,get_drf,rdeqiv, &
	 get_fnet,get_cvf,fnetmaxl

	CONTAINS
!----------------------------------------------------------------
SUBROUTINE get_bbf( temp, bf )	! 黒体放射束密度の計算

	REAL,DIMENSION(:),INTENT(in)	:: temp	! 温度
	REAL,DIMENSION(:),INTENT(out)	:: bf	! 黒体放射束密度

	bf = sig*temp**4

END SUBROUTINE get_bbf
!----------------------------------------------------------------
FUNCTION TF(tau1,tau2)	! 透過関数

	REAL				:: TF
	REAL,INTENT(in)			:: tau1,tau2

	TF = EXP( ( tau1 - tau2 )*3/2 )

END FUNCTION TF
!----------------------------------------------------------------
SUBROUTINE get_urf1( nu, tau, temp, Fu )! 上向き放射束密度の計算

	INTEGER,INTENT(in)		:: nu	! 大気層番号
	REAL,DIMENSION(0:),INTENT(in)	:: tau	! 光学的深さ
	REAL,DIMENSION(0:),INTENT(in)	:: temp	! 温度
	REAL,DIMENSION(:),INTENT(out)	:: Fu	! 上向き放射束密度

	REAL,DIMENSION(:),ALLOCATABLE	:: bf	! 黒体放射束密度
	REAL 				:: dbdtau
	REAL				:: dtau
	INTEGER 			:: nnl

	nnl = SIZE(tau)-2

	ALLOCATE( bf(0:nnl+1) )
	CALL get_bbf( temp, bf )

	dtau = tau(nu-1) - tau(nu+1)
	dbdtau = ( bf(nu-1) - bf(nu+1) )/dtau

	Fu(nu) = bf(1)*TF( tau(nu), tau(1) ) &
	  	+ bf(nu)*( 1.- TF( tau(nu), tau(1) ) ) &
		+ dbdtau * (( tau(nu) - tau(1) ) * TF( tau(nu), tau(1) ) &
		+ 2/3*( 1. - TF( tau(nu), tau(1) ) ))

	DEALLOCATE(bf)

END SUBROUTINE get_urf1
!----------------------------------------------------------------
SUBROUTINE get_urf( tau, temp, Fu)	! 各層の上向き放射束密度の計算

	REAL,DIMENSION(0:),INTENT(in)	:: tau	! 光学的深さ
	REAL,DIMENSION(0:),INTENT(in)	:: temp	! 温度
	REAL,DIMENSION(:),INTENT(out)	:: Fu	! 上向き放射束密度

	REAL,DIMENSION(:),ALLOCATABLE	:: bf	! 黒体放射束密度
	REAL 				:: dbdtau
	REAL 				:: dtau
	INTEGER 			:: i,nnl

	nnl = SIZE(tau)-2

	ALLOCATE( bf(0:nnl+1) )
	CALL get_bbf( temp, bf )

	DO i = 1,nnl
	   dtau = tau(i-1) - tau(i+1)
	   dbdtau = ( bf(i-1) - bf(i+1) )/dtau

	   Fu(i) = bf(1)*TF( tau(i), tau(1) ) &
		+ bf(i)*( 1.- TF( tau(i), tau(1) ) ) &
		+ dbdtau * ( ( tau(i) - tau(1) ) * TF( tau(i), tau(1) ) &
		+ 2/3*( 1. - TF( tau(i), tau(1) ) ) )
	END DO

	DEALLOCATE(bf)

END SUBROUTINE get_urf
!----------------------------------------------------------------
SUBROUTINE get_drf1( nu, tau, temp, Fd )	! 下向き放射束密度の計算

	INTEGER,INTENT(in)		:: nu	! 大気層番号
	REAL,DIMENSION(0:),INTENT(in)	:: tau	! 光学的深さ
	REAL,DIMENSION(0:),INTENT(in)	:: temp	! 温度
	REAL,DIMENSION(:),INTENT(out)	:: Fd	! 下向き放射束密度

	REAL,DIMENSION(:),ALLOCATABLE	:: bf	! 黒体放射束密度
	REAL 				:: dbdtau
	REAL 				:: dtau
	INTEGER 			:: nnl

	nnl = SIZE(tau)-2
	ALLOCATE( bf(0:nnl+1) )
	CALL get_bbf( temp, bf )

	dtau = tau(nu-1) - tau(nu+1)
	dbdtau = ( bf(nu-1) - bf(nu+1) )/dtau

	Fd(nu) = bf(nnl)*TF( tau(nnl), tau(nu) ) &
		+ bf(nu)*( 1. - TF( tau(nnl), tau(nu) ) ) &
		+ dbdtau*( ( tau(nu) - tau(nnl) )*TF( tau(nnl), tau(nu) )  &
		- 2/3*( 1. - TF( 0., tau(nu) ) ) )

	DEALLOCATE(bf)

END SUBROUTINE get_drf1
!----------------------------------------------------------------
SUBROUTINE get_drf(tau,temp,Fd)	! 各層の下向き放射束密度の計算

	REAL,DIMENSION(0:),INTENT(in)	:: tau	! 光学的深さ
	REAL,DIMENSION(0:),INTENT(in)	:: temp	! 温度
	REAL,DIMENSION(:),INTENT(out)	:: Fd	! 下向き放射束密度

	REAL,DIMENSION(:),ALLOCATABLE	:: bf	! 黒体放射束密度
	REAL 				:: dbdtau
	REAL 				:: dtau
	INTEGER 			:: i,nnl

	nnl = SIZE(tau) - 2
	ALLOCATE( bf(0:nnl+1) )
	CALL get_bbf( temp, bf )

	DO i = 1,nnl
	   dtau = tau(i-1) - tau(i+1)
	   dbdtau = ( bf(i-1) - bf(i+1) )/dtau

	   Fd(i) = bf(nnl)*TF( tau(nnl), tau(i) ) &
		+ bf(i)*( 1. - TF( tau(nnl), tau(i) ) ) &
		+ dbdtau*( ( tau(nnl) - tau(i) )*TF( tau(nnl), tau(i) )  &
		- 2/3*( 1. - TF( tau(nnl), tau(i) ) ) )
	END DO

	DEALLOCATE(bf)

END SUBROUTINE get_drf
!----------------------------------------------------------------
SUBROUTINE get_cvf( Fu, Fd, Fsol, Fcnv )	! 対流熱流速密度の計算

	REAL,DIMENSION(:),INTENT(in)	:: Fu	! 上向き射束密度
	REAL,DIMENSION(:),INTENT(in)	:: Fd	! 下向き射束密度
	REAL,DIMENSION(:),INTENT(in)	:: Fsol	! 太陽放射束密度
	REAL,DIMENSION(:),INTENT(out)	:: Fcnv	! 対流熱流束密度

	Fcnv = Fd + Fsol - Fu

END SUBROUTINE get_cvf
!----------------------------------------------------------------
SUBROUTINE get_fnet( Fu, Fd, Fnet )	! 正味放射束密度の計算

	REAL,DIMENSION(:),INTENT(in)	:: Fu	! 上向き射束密度
	REAL,DIMENSION(:),INTENT(in)	:: Fd	! 下向き射束密度
	REAL,DIMENSION(0:),INTENT(out)	:: Fnet	! 正味上向き放射束密度

	Fnet(1:) = Fu(1:) - Fd(1:)
	Fnet(0) = Fnet(1)
	Fnet(SIZE(Fnet)-1) = Fnet(SIZE(Fnet)-2)

END SUBROUTINE get_fnet
!----------------------------------------------------------------
SUBROUTINE rdeqiv( TP, tau, bf, temp, Ftop)	! 放射平衡解の計算

	INTEGER,INTENT(in)		:: TP	! 成層圏下端高度番号
	REAL,DIMENSION(0:),INTENT(in)	:: tau	! 光学的厚さ
	REAL,DIMENSION(0:),INTENT(inout):: bf	! 黒体放射束密度
	REAL,DIMENSION(0:),INTENT(inout):: temp	! 温度
	REAL,INTENT(out) :: Ftop		! 大気上端上向き放射束密度

!-大気上端上向き放射束密度の計算

	Ftop = 2*sig*temp(TP)**4/( 1.5*tau(TP) + 1. )

!-成層圏の温度; 放射平衡解

  	bf(TP+1:)   = 0.5*Ftop*( 1.5*tau(TP+1:) + 1. )
	temp(TP+1:) = ( bf(TP+1:)/sig )**(0.25)

END SUBROUTINE rdeqiv
!----------------------------------------------------------------
SUBROUTINE fnetmaxl( Fnet, LTRP )	! 正味放射束密度の計算

	REAL,DIMENSION(0:),INTENT(in)	:: Fnet	! 正味上向き放射束密度
	INTEGER,INTENT(inout)		:: LTRP	! 圏界面高度番号
	INTEGER				:: nnl	!
	INTEGER				:: i	!

	nnl = SIZE(Fnet)-1

	DO i = nnl,1,-1
	   IF ( Fnet(i) >= Fnet(LTRP) ) THEN
	      LTRP = i
	   END IF
	END DO

END SUBROUTINE fnetmaxl
!----------------------------------------------------------------
END MODULE rad_module