!=================================================================
!  正味放射フラックス日変化の描画
!
!  1999/05/31  小高正嗣
!  1999/06/04  小高正嗣
!  2000/03/20  小高正嗣
!  2000/04/19  小高正嗣
!  2000/05/25  小高正嗣
!=================================================================
      program rad
!-----------------------------------------------------------------      
      implicit real*8 ( a-h,o-z )
!-----------------------------------------------------------------      
#include "grid_size_M.f"
!-----------------------------------------------------------------      
      parameter ( NBND = 16 )
      parameter ( NBNDSR1= 10 ) ! NIR band number 1(4.3μm)
      parameter ( NBNDSR2= 10 ) ! NIR band number 1(2.7μm)
      parameter ( NBNDSR3=  8 ) ! NIR band number 1(2.0μm)
      parameter ( DELNU = 25.0D0 ) ! band width
      parameter ( SIGMA = 5.67D-8 )
      parameter ( PI = 3.14159263D0 )
      parameter ( TSOL   = 5760.0D0 ) ! surface temp. of Sun
      parameter ( DAY    = 86400.0D0)
      parameter ( HOUR   = 3600.0D0 ) 
     

      dimension
     \   ZGRID ( -2:NZB ), TEMP0 ( -2:NZB ), TPOT0 ( -2:NZB ),
     \   DENS0 ( -2:NZB ), PPAI0 ( -2:NZB ), PRES0 ( -2:NZB ),
     \   RVAP0 ( -2:NZB ), QVAP0 ( -2:NZB )

      dimension
     \  TPOT ( -2:NXB, -2:NZB, NROT ), TEMP ( -2:NXB, -2:NZB ), 
     \  QVAP ( -2:NXB, -2:NZB, NROT )

      dimension
     \  tempg ( -2:nxb,nzg )

      DIMENSION
     \   FDZ  ( -2:NZB ), FDZM ( -2:NZB )

      DIMENSION
     \   DTPRAD ( -2:NXB, -2:NZB ),
     \   DTIRAD ( -2:NXB, -2:NZB ),
     \   DTSRAD ( -2:NXB, -2:NZB ),
     \   DTQIRAD( -2:NXB, -2:NZB ),
     \   DTQSRAD( -2:NXB, -2:NZB )

      DIMENSION
     \   OPL   (-2:NZB), OPL_S(-2:NZB)

      dimension 
     \   QTAU   (-2:NXB,-2:NZB), QDTAU   (-2:NXB,-2:NZB),
     \   QTAUIR1(-2:NXB,-2:NZB), QDTAUIR1(-2:NXB,-2:NZB),
     \   QTAUIR2(-2:NXB,-2:NZB), QDTAUIR2(-2:NXB,-2:NZB)

      DIMENSION
     \   TAUN    (0:NZ+1,0:NZ+1,NBND),
     \   TAUN_SR1(0:NZ+1,0:NZ+1,NBNDSR1),
     \   TAUN_SR2(0:NZ+1,0:NZ+1,NBNDSR2),
     \   TAUN_SR3(0:NZ+1,0:NZ+1,NBNDSR3)

      dimension
     \   WNU(NBND), SNU(NBND), RNU(NBND)

      dimension
     \   WNU_SR1(NBNDSR1), SNU_SR1(NBNDSR1), RNU_SR1(NBNDSR1),
     \   WNU_SR2(NBNDSR2), SNU_SR2(NBNDSR2), RNU_SR2(NBNDSR2),
     \   WNU_SR3(NBNDSR3), SNU_SR3(NBNDSR3), RNU_SR3(NBNDSR3)

      dimension
     \   BBRAD ( -2:NXB, -2:NZB ), FXRADU( -2:NXB, -2:NZB ),
     \   FXRADD( -2:NXB, -2:NZB ), FXRADN( -2:NXB, -2:NZB ),
     \   BBRAD_S( -2:NXB )

      dimension
     \   SXRADU (-2:NXB,-2:NZB), SXRADD (-2:NXB,-2:NZB), 
     \   SXRADN (-2:NXB,-2:NZB), 
     \   SXRADU1(-2:NXB,-2:NZB), SXRADD1(-2:NXB,-2:NZB), 
     \   SXRADN1(-2:NXB,-2:NZB), 
     \   SXRADU2(-2:NXB,-2:NZB), SXRADD2(-2:NXB,-2:NZB), 
     \   SXRADN2(-2:NXB,-2:NZB), 
     \   SXRADU3(-2:NXB,-2:NZB), SXRADD3(-2:NXB,-2:NZB), 
     \   SXRADN3(-2:NXB,-2:NZB), 
     \   SQRAD1 (-2:NXB,-2:NZB), 
     \   SQRAD2 (-2:NXB,-2:NZB), 
     \   SQRAD3 (-2:NXB,-2:NZB),
     \   F0(-2:NXB), FXRADG(-2:NXB,-2:NZB), F1(-2:NXB)

      DIMENSION
     \   FXRADU_NU(-2:NXB,-2:NZB,NBND), 
     \   FXRADD_NU(-2:NXB,-2:NZB,NBND),
     \   SXRADU_NU1(-2:NXB,-2:NZB,NBNDSR1), 
     \   SXRADD_NU1(-2:NXB,-2:NZB,NBNDSR1),
     \   SXRADU_NU3(-2:NXB,-2:NZB,NBNDSR3), 
     \   SXRADD_NU3(-2:NXB,-2:NZB,NBNDSR3)

      dimension
     \   FXQRADU (-2:NXB,-2:NZB), FXQRADD (-2:NXB,-2:NZB), 
     \   FXQRADN (-2:NXB,-2:NZB),
     \   FXQRADU1(-2:NXB,-2:NZB), FXQRADD1(-2:NXB,-2:NZB), 
     \   FXQRADN1(-2:NXB,-2:NZB), 
     \   FXQRADU2(-2:NXB,-2:NZB), FXQRADD2(-2:NXB,-2:NZB), 
     \   FXQRADN2(-2:NXB,-2:NZB),
     \   QRAD1   (-2:NXB,-2:NZB), QRAD2   (-2:NXB,-2:NZB)

      dimension
     \   SXQRADU(-2:NXB,-2:NZB), SXQRADD(-2:NXB,-2:NZB), 
     \   SXQRADN(-2:NXB,-2:NZB), SXSOL  (-2:NXB,-2:NZB)

      dimension
     \   SCRADW1(-2:NXB,-2:NZB), SCRADW2(-2:NXB,-2:NZB) 

      dimension 
     \   ASR(0:NZ,0:NZ), BSR(0:NZ,0:NZ), CSR(0:NZ,0:NZ), 
     \   DSR(0:NZ,0:NZ), RSR(0:NZ)     , SSR(0:NZ)

      dimension
     \  FXRADN_mean(-2:NZB)

      dimension TSFC( -2:NXB ), QVSFC( -2:NXB ), FCOLI( -2:NXB )
      
      real*4 F0MEAN, FXRADD_M, FXRADG_M, FXQRADD_M, F1MEAN
      real*4 FXQRADD1_M, FXQRADD2_M

      character*5 header
      character*2 date

      character*100 dumy

      NAMELIST /CONST1/ GRAV, PSD, AMDRY, RDRY, CPDRY, RMVAP, ALATNT
      NAMELIST /CONST2/ CDRAG, VSFC0, PSFC

!-----------------------------------------------------------------      
! 物理定数の設定
!-----------------------------------------------------------------      

      call SETCST
     O   ( GRAV, PSD, AMDRY, RDRY, CPDRY, RMVAP, ALATNT )

      CALL SETCS2
     O   ( CDRAG, VSFC0, TSFC(-2), QVSFC(-2),
     O     FCOLI,
     O     PSFC,
     C     RMVAP,
     D     NX, NXB                      )

      call SETCSP
     C   ( S00, EXCTL, OBLIC, OMGA0 )

      call SETDUST
     O   ( QEXT_SR, QEXT_IR1, QEXT_IR2, 
     O     QALB_SR, QALB_IR1, QALB_IR2, 
     O     QASM_SR, QASM_IR1, QASM_IR2,
     O     QDENS  , QREFF   , QRMOD     )
C
C *** Delta-Eddington Approximation *** (2000/11/09)
C
      QALB_SRE = QALB_DE( QALB_SR, QASM_SR ) 
      QASM_SRE = QASM_DE( QASM_SR )

      QALB_IR1 = QALB_DE( QALB_IR1, QASM_IR1 )
      QASM_IR1 = QASM_DE( QASM_IR1 )

      QALB_IR2 = QALB_DE( QALB_IR2, QASM_IR2 )
      QASM_IR2 = QASM_DE( QASM_IR2 )

!-----------------------------------------------------------------      

      open (11,file='exparam')
      read (11,*) irun
      read (11,*) time0, nloop1, nloop2, dtime, ntmoni, ntkubu, TIMES
      read (11,*) YLAT, ALS
      close(11)

      YLAT = YLAT * PI/180.0D0  ! transfer from degree to rad
      ALS  = ALS  * PI/180.0D0  ! transfer from degree to rad

      tinterv = nloop2 

      call FOPEN ( irun+1 )

      do iz = 1,9
         read(16,'(a100)') dumy
      end do

      do iz = -2, nzb
         read(16,*) ii, zgrid(iz), pres0(iz), dens0(iz), temp0(iz),
     \                     tpot0(iz), ppai0(iz), qvap0(iz), rvap0(iz)
      end do

!-----------------------------------------------------------------      
! 放射場の初期条件設定; バンドパラメターの設定
! (Houghton, 1986;The physics of atmosphere, 2nd ed., 
!  Cambridge Univ. Press)

      DATA ( WNU(I), I=1,NBND )
     1     /  512.5D0,  537.5D0,  562.5D0,  587.5D0,
     2        612.5D0,  637.5D0,  662.5D0,  687.5D0,  
     3        712.5D0,  737.5D0,  762.5D0,  787.5D0,
     4        812.5D0,  837.5D0,  862.5D0,  887.5D0  / 

      DATA ( SNU(I), I=1,NBND )
     1     /  1.952D-2,  2.785D-1,  5.495D-1,  5.331D0, 
     2        5.196D+2,  7.778D+3,  8.746D+4,  2.600D+4,
     3        1.232D+3,  2.042D+2,  7.278D0 ,  1.337D0,
     4        3.974D-1,  1.280D-2,  2.501D-3,  3.937D-3  /

      DATA ( RNU(I), I=1,NBND )
     1     /  2.870D-1,  1.215D0 ,  2.404D0 ,  1.958D+1, 
     2        5.804D+1,  2.084D+2,  7.594D+2,  2.635D+2,
     3        8.387D+1,  2.852D+1,  6.239D0 ,  2.765D0 ,
     4        8.897D-1,  3.198D-1,  1.506D-1,  1.446D-1  /

! CO2 4.3 μm

      DATA ( WNU_SR1(I), I=1,NBNDSR1 )
     1     /  2212.5D0, 2237.5D0, 2262.5D0, 2287.5D0,
     2        2312.5D0, 2337.5D0, 2362.5D0, 2387.5D0,  
     3        2412.5D0, 2437.5D0                         /

      DATA ( SNU_SR1(I), I=1,NBNDSR1 )
     1     /  9.504D-1,  2.217D+2,  4.566D+3,  7.965D+3, 
     2        1.055D+5,  5.587D+5,  6.819D+5,  1.256D+4,
     3        7.065D-2,  8.522D-2                        /

      DATA ( RNU_SR1(I), I=1,NBNDSR1 )
     1     /  2.866D0 ,  3.000D+1,  1.134D+2,  2.011D+2, 
     2        5.880D+2,  1.206D+3,  1.182D+3,  8.873D+1,
     3        3.404D-1,  4.236D-1                        /

! CO2 2.7 μm

      DATA ( WNU_SR2(I), I=1,NBNDSR2 )
     1     /  3150.0D0, 3250.0D0, 3350.0D0, 3450.0D0,
     2        3550.0D0, 3650.0D0, 3750.0D0, 3850.0D0,  
     3        3950.0D0, 4050.0D0                         /

      DATA ( SNU_SR2(I), I=1,NBNDSR2 )
     1     /  1.324D-1,  7.731D-2,  1.232D0 ,  5.159D0, 
     2        4.299D+3,  1.543D+4,  1.649D+4,  1.180D-1,
     3        1.464D-2,  1.251D-2                        /

      DATA ( RNU_SR2(I), I=1,NBNDSR2 )
     1     /  9.836D-1,  4.900D-1,  2.952D0 ,  7.639D0 , 
     2        1.914D+2,  3.245D+2,  2.722D+2,  9.535D-1,
     3        2.601D-1,  2.021D-1                        /


! CO2 2.0 μm

      DATA ( WNU_SR3(I), I=1,NBNDSR3 )
     1     /  4650.0D0, 4750.0D0, 4850.0D0, 4950.0D0,
     2        5050.0D0, 5150.0D0, 5250.0D0, 5350.0D0     /

      DATA ( SNU_SR3(I), I=1,NBNDSR3 )
     1     /  2.185D-1,  2.040D0 ,  1.197D+2,  4.829D+2, 
     2        8.778D+1,  8.346D+1,  8.518D-2,  4.951D-1  /

      DATA ( RNU_SR3(I), I=1,NBNDSR3 )
     1     /  1.916D0,   6.475D0 ,  3.112D+1,  5.759D+1, 
     2        2.012D+1,  1.804D+1,  8.474D-1,  1.597D0   / 

!-----------------------------------------------------------------      

      CALL SETGRD
     O   ( ZGRID(-2), FDZ(-2), FDZM(-2), DX,
     D     NZ, NZB                           )

      call RBASIC
     I   ( PRES0, PSFC , SNU  , RNU  , 
     O     OPL  , TAUN ,
     C     GRAV , DELNU,
     D     NZ   , NZB  , NBND  )

!-----------------------------------------------------------------      

      write(6,*) 'data number = ', nloop1

      do it = 1, nloop1

         read(34) TIME, ((TPOT(IX,IZ,1),IX=-2,nxb), IZ=-2,NZB )
         read(35) time, ((qvap(ix,iz, 1),ix=-2,nxb), iz=1,nzb)
         read(40) TIME, ((tempg(ix,iz),ix=-2,nxb), iz=1,nzg)
         write(0,*) 'reading time = ', time
!-----------------------------------------------------------------      

         do IX = -2,NXB
            TSFC(IX) = TEMPG(ix,NZG)
         end do

         call CLFXSOL
     I      ( S00  , ALS   , YLAT , TIME, TIMES,
     O        S0   , COSZ  , 
     C        EXCTL, OBLIC , OMGA0, DAY , PI    )

C
C *** CO2 Infrared Heaating/Cooling ***
C
         CALL CLRAD
     I      ( TPOT  , TAUN  , PSFC  ,
     I        1,
     O        DTIRAD, 
     B        TPOT0 , TEMP0 , PRES0 , DENS0 , FDZ  ,
     C        WNU   ,
     C        TSFC  , CPDRY , RDRY  , DX    , DELNU,  GRAV,
     W        TEMP  , BBRAD , BBRAD_S, 
     W        FXRADU, FXRADD, FXRADN, FXRADU_NU, FXRADD_NU,
     D        NX    , NXB   , NZ    , NZB   , NROT  , NBND  )

C
C *** Dust Optical Depth Set up *** (2000/11/09)
C

         do ITQ = 1, NROT
         do IZ = -2, NZB
         do IX = -2, NXB
            QVAP(IX,IZ,ITQ) = QVAP0(IZ)
         end do
         end do
         end do

         call SETQTAU
     I      ( DENS0, QVAP , 
     I        1, 
     O        QTAU , QDTAU,
     O        QTAUIR1, QDTAUIR1, QTAUIR2, QDTAUIR2,
     C        QREFF  , QDENS   , QALB_SR, QASM_SR , 
     C        QEXT_SR, QEXT_IR1, QEXT_IR2, PI     ,
     G        FDZ  , DX,
     D        NX   , NXB  , NZ   , NZB   , NROT     )

C         write(*,*) QTAU(0,0), QTAUIR1(0,0), QTAUIR2(0,0)

C
C *** Dust Infrared Radiative Heating/Cooling *** (2000/11/09)
C For 5-11.6 micro m 
C
         DINU1  = 900.0D0
         DINU2  = 2000.0D0
         DELNUD = DELNU * 4.0 

         call CLDIRAD
     I     (  TPOT  , PSFC  , 1,
     I        QTAUIR1, QDTAUIR1, 
     I        DINU1  , DINU2   , DELNUD,
     O        QRAD1  , 
     B        PRES0  , DENS0   , TEMP0, TPOT0, 
     C        QALB_IR1, QASM_IR1, 
     C        RDRY  , CPDRY , TSFC  , GRAV  , PI   , 
     W        TEMP  , BBRAD , BBRAD_S, 
     W        FXQRADU1, FXQRADD1, FXQRADN1, 
     W        SCRADW1, SCRADW2, 
     W        ASR   , BSR   , CSR   , DSR  ,
     W        RSR   , SSR   ,
     G        FDZ   , FDZM  , DX    ,
     D        NX    , NXB   , NZ    , NZB  , NROT      )

C
C For 20-200 micro m 
C
         DINU1  = 50.0D0
         DINU2  = 500.0D0
         DELNUD = DELNU * 2.0 

         call CLDIRAD
     I     (  TPOT  , PSFC  , 1,
     I        QTAUIR2, QDTAUIR2, 
     I        DINU1  , DINU2   , DELNUD,
     O        QRAD2  , 
     B        PRES0  , DENS0   , TEMP0, TPOT0, 
     C        QALB_IR2, QASM_IR2, 
     C        RDRY  , CPDRY , TSFC  , GRAV  , PI   , 
     W        TEMP  , BBRAD , BBRAD_S, 
     W        FXQRADU2, FXQRADD2, FXQRADN2, 
     W        SCRADW1, SCRADW2, 
     W        ASR   , BSR   , CSR   , DSR  ,
     W        RSR   , SSR   ,
     G        FDZ   , FDZM  , DX    ,
     D        NX    , NXB   , NZ    , NZB  , NROT      )

         do iz = -2, NZB
         do ix = -2, NXB
C            DTQIRAD(ix,iz) = QRAD1   (ix,iz) + QRAD2   (ix,iz)
C            FXQRADU(ix,iz) = FXQRADU1(ix,iz) + FXQRADU2(ix,iz)
            FXQRADD(ix,iz) = FXQRADD1(ix,iz) + FXQRADD2(ix,iz)
C            FXQRADN(ix,iz) = FXQRADN1(ix,iz) + FXQRADN2(ix,iz)
         end do
         end do

C
C *** Solar Heating/Cooling *** (2000/10/31)
C
         if ( COSZ .gt. 0.0D0 ) then

            do IZ = -2, NZB
               OPL_S(IZ) = OPL(IZ) / COSZ
            end do
         
            call CLTAUN
     I         ( OPL_S, SNU_SR1, RNU_SR1, PRES0,
     O           TAUN_SR1, 
     C           DELNU, 
     D           NZ   , NZB, NBNDSR1 )

            call CLTAUN
     I         ( OPL_S, SNU_SR2, RNU_SR2, PRES0,
     O           TAUN_SR2, 
     C           DELNU*4.0, 
     D           NZ   , NZB, NBNDSR2 )

            call CLTAUN
     I         ( OPL_S, SNU_SR3, RNU_SR3, PRES0,
     O           TAUN_SR3, 
     C           DELNU*4.0, 
     D           NZ   , NZB, NBNDSR3 )

C *** CO2 4.3μm ***

            call CLSRAD
     I         ( TAUN_SR1, S0    , COSZ  ,
     O           SXRADU1, SXRADD1, SXRADN1, SQRAD1,
     B           TPOT0 , TEMP0 , PRES0 , DENS0 , FDZ  ,
     C           WNU_SR1, 
     C           CPDRY , TSOL  , DX    , SIGMA, 
     C           DELNU , GRAV  , PI    ,
     W           SXRADU_NU1, SXRADD_NU1,
     D           NX ,  NXB,  NZ  , NZB   , NROT, NBNDSR1 )

C *** CO2 2.7μm ***

            call CLSRAD
     I        ( TAUN_SR2, S0    , COSZ  ,
     O           SXRADU2, SXRADD2, SXRADN2, SQRAD2,
     B           TPOT0 , TEMP0 , PRES0 , DENS0 , FDZ  ,
     C           WNU_SR2, 
     C           CPDRY , TSOL  , DX    , SIGMA, 
     C           DELNU*4.0 , GRAV  , PI    ,
     W           SXRADU_NU1, SXRADD_NU1,
     D           NX ,  NXB,  NZ  , NZB   , NBNDSR2 )

C *** CO2 2.0μm ***

            call CLSRAD
     I         ( TAUN_SR3, S0    , COSZ  ,
     O           SXRADU3, SXRADD3, SXRADN3, SQRAD3,
     B           TPOT0 , TEMP0 , PRES0 , DENS0 , FDZ  ,
     C           WNU_SR3, 
     C           CPDRY , TSOL  , DX    , SIGMA, 
     C           DELNU*4.0 , GRAV  , PI    ,
     W           SXRADU_NU3, SXRADD_NU3,
     D           NX ,  NXB,  NZ  , NZB   , NBNDSR3 )
         
            do iz = -2, NZB
               do ix = -2, NXB
                  DTSRAD(iz,iz) = SQRAD1(ix,iz)  
     \                 + SQRAD2(ix,iz)  + SQRAD3(ix,iz)
                  SXRADD(iz,iz) = SXRADD1(ix,iz) 
     \                 + SXRADD2(ix,iz) + SXRADD3(ix,iz)
               end do
            end do

C *** Dust *** (2000/11/09)

            call CLDSRAD
     I         ( QTAU  , QDTAU  , S0    , COSZ  , 
     O           DTQSRAD, SXSOL , 
     B           TPOT0 , TEMP0  , PRES0 , DENS0 ,
     C           ALBEDO, QALB_SRE, QASM_SRE,
     C           DELNU , CPDRY  , TSOL   , SIGMA, 
     C           GRAV  , PI     , 
     W           SXQRADU, SXQRADD, SXQRADN, 
     W           SCRADW1, SCRADW2,
     W           ASR   , BSR    , CSR    , DSR  ,
     W           RSR   , SSR    ,
     G           FDZ   , FDZM   , DX    ,
     D           NX    , NXB    , NZ     , NZB           )

         else 

            do iz = -2, NZB
               do ix = -2, NXB
                  SXRADU(ix,iz) = 0.0D0
                  SXRADD(ix,iz) = 0.0D0
                  SXRADN(ix,iz) = 0.0D0
                  SXQRADU(ix,iz)= 0.0D0
                  SXQRADD(ix,iz)= 0.0D0
                  SXQRADN(ix,iz)= 0.0D0
                  SXSOL (ix,iz) = 0.0D0
                  SQRAD1(ix,iz) = 0.0D0
                  SQRAD2(ix,iz) = 0.0D0 
                  SQRAD3(ix,iz) = 0.0D0
                  DTSRAD(ix,iz) = 0.0D0
                  DTQSRAD(ix,iz)= 0.0D0
               end do
            end do

         end if

         do ix = -2, NXB
            F0(ix) = ( S0 - SXSOL(ix,NZ+1) ) * COSZ + SXQRADD(ix,0) - 
     \           ( SXRADD (ix,NZ+1) - SXRADD(ix,0)  ) 
         end do

         do iz = -2, NZB
            do ix = -2, NXB
               FXRADG(ix,iz) = FXRADD(ix,iz) + FXQRADD(ix,iz)
            end do
         end do
         
         F0MEAN = 0.0
         F1MEAN = 0.0
         FXRADD_M = 0.0
         FXRADG_M = 0.0
         FXQRADD_M = 0.0
         FXQRADD1_M = 0.0
         FXQRADD2_M = 0.0

         do ix = 1, NX
            F0MEAN = F0MEAN + F0(IX) / NX
            F1MEAN = F1MEAN + F1(IX) / NX
            FXRADG_M = FXRADG_M + FXRADG(IX,0) / NX
            FXRADD_M = FXRADD_M + FXRADD(IX,0) / NX
            FXQRADD_M= FXQRADD_M + FXQRADD(IX,0) / NX
            FXQRADD1_M= FXQRADD1_M + FXQRADD1(IX,0) / NX
            FXQRADD2_M= FXQRADD2_M + FXQRADD2(IX,0) / NX
         end do

         write(*,*) FXRADG_M, FXRADD_M, FXQRADD1_M, FXQRADD2_M
C         write(*,*) F0MEAN, F1MEAN
            
      end do


      stop
      end

