C=================================================================
C
C     PROGRAM netheat
C
C     正味加熱量の日変化図作成プログラム
C
C     履歴 2000/05/25 小高正嗣
C          2000/11/21 小高正嗣; D 論用に書き直し
C=================================================================
      program netheat
!-----------------------------------------------------------------      
      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 ) 

      parameter ( NT = 100 )
     
      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 ), CDTURB( -2:NXB, -2:NZB, NROT ),
     \  CDZ  ( -2:NZB )

      dimension
     \     FUMOMA(-2:NXB,2), FVMOMA(-2:NXB,2), FHEATA(-2:NXB,2), 
     \     FEVAPA(-2:NXB,2), FWPREA(-2:NXB,2)

      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 ),
     \   DTPVIS ( -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)

      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 TSFC( -2:NXB ), QVSFC( -2:NXB ), FCOLI( -2:NXB )
      
      real*4 
     \     NETIR_G(0:NT), NETSR_G(0:NT), NETIR_A(0:NT), NETSR_A(0:NT),
     \     FHEAT(0:NT), TIMEL(0:NT), FHEAT_G(0:NT), NETVIS_A(0:NT),
     \     HEAT_A(0:NT),HEAT_G(0:NT)

      real*4
     \     YMIN1,  YMAX1,  YMIN2,  YMAX2

      character*14 dheader
      character*2 date

      character*100 dumy

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

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

      YMAX1 = 600.0
      YMIN1 = -YMAX1

      YMAX2 = 50
      YMIN2 = -YMAX2

      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 * ( dtime / 2.)

c*** 時刻変換: 日付設定 ************************************
      time = time0 + times
      write(date, '(I2.2)') INT( time/day ) + 1
      dheader = 'day=' // date 
c**************************************************************

      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

      ilate = INT( TIMES / HOUR )

      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=-2,nzb)
         read(38) TIME, ((cdturb(ix,iz,1),ix=-2,nxb), IZ=-2,nzb)
         read(40) TIME, ((tempg(ix,iz), ix=-2,nxb), IZ=1,nzg)

         read(72) TIME, FUMOMA, FVMOMA, FHEATA, FEVAPA, FWPREA

         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 Heating/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

C         do ITQ = 1, NROT
C         do IZ = -2, NZB
C         do IX = -2, NXB
C            QVAP(IX,IZ,ITQ) = QVAP0(IZ)
C         end do
C         end do
C         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)
C            FXQRADD(ix,iz) = FXQRADD1(ix,iz) + FXQRADD2(ix,iz)
            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(ix,iz) = SQRAD1(ix,iz)  
     \                 + SQRAD2(ix,iz)  + SQRAD3(ix,iz)
                  SXRADD(ix,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
C====================================================================
C Turbulent Disspative Heating

         CMOM=0.2D0
         CDIS=0.2D0

         DO IZ = 0, NZ
            DO IX = 0, NX
               DTPVIS(IX,IZ) = DENS0(iz) * FDZ(iz) * DX * 
     \              CDIS * CDTURB(IX,IZ,1)**3 / CMOM**3 
     \              / ( DX * FDZ(IZ) )**4 
            END DO
         END DO

C
         CALL BOUND ( DTPVIS, NX, NZ, NXB, NZB, 1, 1, 1, 1 )

C====================================================================

         do ix = 1, NX
            FHEAT(it+ilate) = FHEAT(it+ilate) + 
     \           FHEATA(IX, 1) / tinterv / NX 
         end do

         do IX = 1, NX
            NETIR_G(it+ilate) = NETIR_G(it+ilate) + 
     \        REAL( - FXRADN(IX,0) + FXQRADD(IX,0) )/NX

            NETSR_G(it+ilate) = NETSR_G(it+ilate) + 
     \        ( 1.0 - ALBEDO ) * 
     \        REAL( ( S0 - SXSOL(IX,NZ+1) ) * COSZ + SXQRADD(IX,0) - 
     \        ( SXRADD (IX,NZ+1) - SXRADD(IX,0) ) ) / NX 

            NETIR_A(it+ilate) = NETIR_A(it+ilate) +  
     \           ( REAL( ( FXRADN(IX,0) - FXRADN(IX,NZ+1) ) 
     \           + ( FXQRADN(IX,0) - FXQRADN(IX,NZ+1) ) ) ) / NX 

            NETSR_A(it+ilate) = NETSR_A(it+ilate) +
     \           ( REAL( ( SXRADD(IX,NZ+1) - SXRADD(IX,0) ) ) 
     \           + ( SXQRADN(IX,0) - SXQRADN(IX,NZ+1) ) ) / NX

         end do

         do iz = 0, NZ
            CDZ ( iz ) = 0.0D0
            do ix = 1, NX
               CDZ( iz ) = CDZ( iz ) + CDTURB( ix, iz, 1 )**3.0 / NX
            end do
         end do

         NETVIS_A(it+ilate) = 0.0
         do IZ = 0, NZ-1
            NETVIS_A(it+ilate) = NETVIS_A(it+ilate) + 
     \           DENS0(iz) * FDZ(iz) * DX * 
     \           CDIS * CDZ(IZ) / CMOM**3  / ( DX * FDZ(IZ) )**4 
         end do


         TIMEL(it+ilate) = (MOD ( TIME, DAY ) + TIMES ) / HOUR

         write(*,*) TIMEL(it+ilate),
C     \        NETIR_G(it+ilate), NETSR_G(it+ilate), 
C     \        NETIR_A(it+ilate), NETSR_A(it+ilate),
C     \        FHEAT  (it+ilate)
     \         SXRADD(0,NZ+1), SXRADD(0,0), NETVIS_A(it+ilate)

      end do

      NETIR_G(ilate) = NETIR_G(nloop1+ilate)
      NETSR_G(ilate) = NETSR_G(nloop1+ilate)
      NETIR_A(ilate) = NETIR_A(nloop1+ilate)
      NETSR_A(ilate) = NETSR_A(nloop1+ilate)
      FHEAT(ilate)  = FHEAT(nloop1+ilate)
      TIMEL(ilate)  = TIMEL(nloop1+ilate)
      
      NETVIS_A(ilate) = NETVIS_A(nloop1+ilate)

      il = 0
      do it = 1, nloop1+ilate
         if ( TIMEL(it) .ge. 24.0 ) then
            TIMEL(il)  = TIMEL(it) - 24.0
            NETIR_G(il) = NETIR_G(it)
            NETSR_G(il) = NETSR_G(it)
            NETIR_A(il) = NETIR_A(it)
            NETSR_A(il) = NETSR_A(it)
            NETVIS_A(il) = NETVIS_A(it)
            FHEAT(il) = FHEAT(it)
            il = il + 1
         end if
      end do

      do it = 0, nloop1
         FHEAT_G(it)= - FHEAT(it)
      end do

      do it = 0, nloop1
         HEAT_G(it) = NETIR_G(it) + NETSR_G(it) + FHEAT_G(it)
         HEAT_A(it) = NETIR_A(it) + FHEAT(it)   + NETVIS_A(it)
      end do

      call sgpwsn
      read(*,*) iws
      CALL SWISET('IWIDTH', 840)
      CALL SWISET('IHEIGHT', 420)
      call sglset( 'LFULL', .TRUE. )
      call swlset( 'LDUMP', .TRUE. )
      call gropn(iws)
      call sldiv('Y', 2, 1 )
      call sglset( 'LCORNER', .FALSE. )

      call grfrm
      CALL sglset( 'LCHAR', .FALSE. )
      call grswnd( 0.0, 24.0, YMIN1, YMAX1 )   
      call grsvpt( 0.2, 0.8, 0.2, 0.8)                     
      call grstrn( 1 )
      call grstrf
      call ussttl( 'Local Time', 'Hour', 'Heat Flux', 'W/m^2' )
      call usdaxs

      call SGSPLI( 43 )
      call sgplu(nloop1+1, TIMEL, NETSR_G)
      call SGSPLI( 23 )
      call sgplu(nloop1+1, TIMEL, NETIR_G)
      call SGSPLI( 33 )
      call sgplu(nloop1+1, TIMEL, FHEAT_G)
      call SGSPLI( 773 )
      call sgplu(nloop1+1, TIMEL, HEAT_G)
      call uxsttl( 'T', dheader, 1. )      
     
      call grfrm
      CALL sglset( 'LCHAR', .FALSE. )
      call grswnd( 0.0, 24.0, YMIN2, YMAX2 )   
      call grsvpt( 0.2, 0.8, 0.2, 0.8)                     
      call grstrn( 1 )
      call grstrf
      call ussttl( 'Local Time', 'Hour', 'Heat Flux', 'W/m^2' )
      call usdaxs

      call SGSPLI( 43 )
      call sgplu(nloop1+1, TIMEL, NETSR_A)
      call SGSPLI( 23 )
      call sgplu(nloop1+1, TIMEL, NETIR_A)
      call sgspli( 33 )
      call sgplu(nloop1+1, TIMEL, FHEAT )
      call sgspli( 63 )
      call sgplu(nloop1+1, TIMEL, NETVIS_A)
      call SGSPLI( 773 )
      call sgplu(nloop1+1, TIMEL, HEAT_A)
      call uxsttl( 'T', dheader, 1. )      
     
      call grcls                                                        

      call FCLOSE





      stop
      end

