C********************************************************************
C   PROGRAMME FOR CALCULATING TMRT
C     coded by Pek ( ------- ) 
C     modified by M.Kanda ( 2000 7.21 ) Kimura and Liu equation
C　　 modified by M.Kanda ( 2000 7.24 ) Tmrt 
C********************************************************************
                   IMPLICIT REAL*8(A-H,O-Z)
C
      CALL CONSTANT 
      CALL SUNRAY
      CALL RADIATION
      CALL OUTDAT
C
      STOP 
      END 
C********************************************************************
      SUBROUTINE OUTDAT
C********************************************************************
                   IMPLICIT REAL*8(A-H,O-Z)
	PARAMETER (IMAX=6724,NMAX=25)
      COMMON /RADCOMP/ RSDN(NMAX,IMAX),RSUP(NMAX,IMAX)
     &                ,RLDN(NMAX,IMAX),RLUP(NMAX,IMAX)
      COMMON /RADCST1/ RADIO,ETA,SIGMA,HALB
      COMMON /RADCST2/ TM(NMAX),BETAH(NMAX),AH(NMAX)
      COMMON /RADCST3/ EIDO,EKEIDO,ETOSHI,DAYJUL
      COMMON /RADHUM / HRDIR(NMAX,IMAX),HRDIF(NMAX,IMAX)
     &                ,HLONG(NMAX,IMAX),HRTOTAL(NMAX,IMAX)
     &                ,TMRT(NMAX,IMAX),RDIR(NMAX,IMAX),RDIF(NMAX,IMAX)
      COMMON /PHYSE  / UB(NMAX,IMAX),TB(NMAX,IMAX)
     &                ,RB(NMAX,IMAX),RBH(NMAX,IMAX)
C 
      OPEN(UNIT=30,FILE='TMRT-1',STATUS='UNKNOWN',FORM='FORMATTED')
C
	                  DO 200 N = 1,NMAX
			          DO 200 I = 1,IMAX 
      WRITE(30,123) UB(N,I),TB(N,I),RB(N,I),RBH(N,I),TMRT(N,I)
  200           CONTINUE
  123 FORMAT(5F12.6)
C
	CLOSE(30)
      RETURN
      END
C********************************************************************
       SUBROUTINE CONSTANT
C********************************************************************
                   IMPLICIT REAL*8(A-H,O-Z)
	PARAMETER (IMAX=6724,NMAX=25)
      COMMON /RADCOMP/ RSDN(NMAX,IMAX),RSUP(NMAX,IMAX)
     &                ,RLDN(NMAX,IMAX),RLUP(NMAX,IMAX)
      COMMON /RADCST1/ RADIO,ETA,SIGMA,HALB
      COMMON /RADCST2/ TM(NMAX),BETAH(NMAX),AH(NMAX)
      COMMON /RADCST3/ EIDO,EKEIDO,ETOSHI,DAYJUL
      COMMON /RADHUM / HRDIR(NMAX,IMAX),HRDIF(NMAX,IMAX)
     &                ,HLONG(NMAX,IMAX),HRTOTAL(NMAX,IMAX)
     &                ,TMRT(NMAX,IMAX),RDIR(NMAX,IMAX),RDIF(NMAX,IMAX)
      COMMON /PHYSE  / UB(NMAX,IMAX),TB(NMAX,IMAX)
     &                ,RB(NMAX,IMAX),RBH(NMAX,IMAX)
C---------------------------------
C ETA   射出率
C RADIO 太陽定数
C SIGMA ステファンボルツマン定数
C----------------------------------
C EIDO   緯度
C EKEIDO 経度
C TM     標準時
C ETOSHI 西暦
C DAYJUL ジュリアン日
C 
      EIDO  = 35.0
      EKEIDO= 135.0
      ETOSHI= 1998
      DAYJUL= 220
C
      OPEN(UNIT=10,FILE='GROUND',STATUS='OLD',FORM='FORMATTED')
      RADIO = 1365
      ETA   = 1
      SIGMA = 0.0000000567
      HALB  = 0.3
      PAI   = 4.0*ATAN(1.0)
C
	          DO 100 N = 1,NMAX
	TM(N) =  1.0*(N-1)
  100           CONTINUE
C
      READ(10,*)
	          DO 200 N = 1,NMAX
			  DO 200 I = 1,IMAX 
      READ(10,*) UB(N,I),TB(N,I),RB(N,I),RBH(N,I)
     &           ,RSDN(N,I),RSUP(N,I),RLDN(N,I),RLUP(N,I)
	UB(N,I)=MAX(UB(N,I),0.5)
C      CSANG = 2*PAI/24.0*(TM(N)-6.0)
C      CCS1  = 0.5*PAI
C      IFLG2 = 1 
C      IF (CSANG.GT.CCS1) IFLG2 = 0
C      SANG  = ABS(0.5*PAI-CSANG) 
C      CC1 = 1000*SIN(CSANG)
C      CC2 = 0.0
C      CC3 = 0.5*PAI
C      SANG  = MIN(CC3,SANG)
C      RSDN(N,I) = MAX(CC2,CC1)
C	RSUP(N,I) = 0.1*RSDN(N,I)
C	RLDN(N,I) = 400
C	RLUP(N,I) = 450 
  200           CONTINUE
C
      CLOSE(10)
      RETURN
      END
C********************************************************************
      SUBROUTINE SUNRAY
C********************************************************************
                   IMPLICIT REAL*8(A-H,O-Z)
	PARAMETER (IMAX=6724,NMAX=25)
      COMMON /RADCOMP/ RSDN(NMAX,IMAX),RSUP(NMAX,IMAX)
     &                ,RLDN(NMAX,IMAX),RLUP(NMAX,IMAX)
      COMMON /RADCST1/ RADIO,ETA,SIGMA,HALB
      COMMON /RADCST2/ TM(NMAX),BETAH(NMAX),AH(NMAX)
      COMMON /RADCST3/ EIDO,EKEIDO,ETOSHI,DAYJUL
      COMMON /RADHUM / HRDIR(NMAX,IMAX),HRDIF(NMAX,IMAX)
     &                ,HLONG(NMAX,IMAX),HRTOTAL(NMAX,IMAX)
     &                ,TMRT(NMAX,IMAX),RDIR(NMAX,IMAX),RDIF(NMAX,IMAX)
      COMMON /PHYSE  / UB(NMAX,IMAX),TB(NMAX,IMAX)
     &                ,RB(NMAX,IMAX),RBH(NMAX,IMAX)
C---------------------------------
C
                   DO 1000 N=1,NMAX
                   DO 1000 I=1,IMAX
C ------------------------------------------------------------
C TIME PARAMETERS---------------------------------------------
C ------------------------------------------------------------
C H 　　太陽高度
C RSDN  全天日射量
C RADKT 正規化された全天日射
C RADKD 正規化された直達日射
C EIDO   緯度
C EKEIDO 経度
C TM     標準時
C ETOSHI 西暦
C DAYJUL ジュリアン日
C 
      PAI      = 4.0*ATAN(1.0)
      EEN      = ETOSHI-1968
      EED      = 3.7+0.2596*EEN-1.0*INT((EEN+3)/4)
      EEM      = 0.9856*(DAYJUL-EED)
      EEE      = 12.39+0.0172*(EEN+EEM/360)
      EEV      = EEM + 1.914*SIN(EEM/180*PAI)+0.02*SIN(2*EEM/180*PAI)
      EET1     = EEM - EEV 
      EET2     = ATAN(0.043*SIN(2*(EEV+EEE)/180*PAI)
     &             /(1-0.043*COS(2*(EEV+EEE)/180*PAI)))
      EET      = EET1 + EET2
      FAI      = EIDO/180*PAI
      DELTA    = ASIN(COS((EEV+EEE)/180*PAI)*SIN(-23.4393/180*PAI))
      TT       = (15*(TM(N)-12)+(EKEIDO-135.0)+EET)/180*PAI
      SINH     = SIN(FAI)*SIN(DELTA)+COS(FAI)*COS(DELTA)*COS(TT)
	BETAH(N) = ASIN(SINH) 
      RADKT    = RSDN(N,I)/(RADIO*SINH) 
      RADKC    = 0.5163+0.333*SINH+0.00803*(SINH*SINH)
      IF ((SINH.LE.0.0).OR.(RSDN(N,I).LE.0.0)) GOTO 222
C
      IF(RADKT.GE.RADKC) 
     &      RADKD=-0.43+1.43*RADKT
      IF(RADKT.LT.RADKC) 
     &      RADKD=(2.277-1.258*SINH+0.2396*SINH*SINH)*RADKT**3
C Kimura Equation
      RDIR1= RADKD*RADIO*SINH        
	RDIF1= RSDN(N,I)-RDIR1              
C Liu and Jourdan Equation
      RDIR2=(RSDN(N,I)-0.2710*RADIO*SINH)/(1-0.2939)
      RDIF2=RSDN(N,I)-RDIR2 
C 散乱日射 RDIF
C 直達日射 RDIR
      IF ((RADKT.LT.0.4).AND.(RDIR2.LT.RDIR1)) THEN
	       RDIF(N,I)=RDIF1
		   RDIR(N,I)=RDIR1
		                                     ELSE
	       RDIF(N,I)=RDIF2
		   RDIR(N,I)=RDIR2 
	                                         END IF
	       GOTO 223 
  222        CONTINUE
             RDIF(N,I)  = 0
	       RDIR(N,I)  = 0
	       BETAH(N) = 0
  223        CONTINUE
 1000        CONTINUE
      RETURN
      END 
C********************************************************************
	 SUBROUTINE RADIATION
C********************************************************************
                   IMPLICIT REAL*8(A-H,O-Z)
	PARAMETER (IMAX=6724,NMAX=25)
      COMMON /RADCOMP/ RSDN(NMAX,IMAX),RSUP(NMAX,IMAX)
     &                ,RLDN(NMAX,IMAX),RLUP(NMAX,IMAX)
      COMMON /RADCST1/ RADIO,ETA,SIGMA,HALB
      COMMON /RADCST2/ TM(NMAX),BETAH(NMAX),AH(NMAX)
      COMMON /RADCST3/ EIDO,EKEIDO,ETOSHI,DAYJUL
      COMMON /RADHUM / HRDIR(NMAX,IMAX),HRDIF(NMAX,IMAX)
     &                ,HLONG(NMAX,IMAX),HRTOTAL(NMAX,IMAX)
     &                ,TMRT(NMAX,IMAX),RDIR(NMAX,IMAX),RDIF(NMAX,IMAX)
      COMMON /PHYSE  / UB(NMAX,IMAX),TB(NMAX,IMAX)
     &                ,RB(NMAX,IMAX),RBH(NMAX,IMAX)
C
                   DO 1000 N=1,NMAX
                   DO 1000 I=1,IMAX
C  直達日射量
      IF (BETAH(N).LE.0) THEN 
	AH(N) = 0.0
	                ELSE 
      AH(N)=(2*14.1*1/TAN(BETAH(N))/3.141593+1)/(2*14.1+2)
	                END IF
      HRDIR(N,I)=AH(N)*RDIR(N,I)
C
C　散乱日射量
      HRDIF(N,I)=(RDIF(N,I)+HALB*RSDN(N,I))*0.5
C
C  長波放射量
      HLONG(N,I)=(RLDN(N,I)+RLUP(N,I))*0.5
C
C  総放射量
      HRTOTAL(N,I)=HRDIR(N,I)+HRDIF(N,I)+HLONG(N,I)
C
C　平均放射温度
      TMRT(N,I)=(HRTOTAL(N,I)/(ETA*SIGMA))**0.25-273.15
C
 1000                 CONTINUE
      RETURN
      END
