C**********************************************************************
C
      PROGRAM DUSTYMARS
C
C     履歴  1999/02/02  小高正嗣; バルク係数パラメタリゼーション
C           1999/06/07  小高正嗣; 放射スキーム組み込み
C           1999/06/14  小高正嗣; 出力データ変更
C           1999/07/21  小高正嗣; 地表面温度計算を放射スキームに組み込み
C           1999/07/23  小高正嗣; 地表面温度計算を指定
C           2000/03/30  小高正嗣; 放射スキームを改良
C           2000/04/18  小高正嗣; 放射を time split で計算
C           2000/05/12  小高正嗣; 射出率スキームへ
C           2000/05/16  小高正嗣; バンドモデルへ(CLRAD)
C           2000/05/25  小高正嗣; コードを整理
C           2000/10/31  小高正嗣; ・温位の非線形数値拡散(NLDIFV)を解除
C                                 ・地表面温度を陽に計算(CLTSFC)
C                                 ・近赤外放射を計算(CLSRAD)
C           2000/11/09  小高正嗣; ・ダストの放射計算(CLQIRAD, CLQSRAD)
C                                 ・ダスト分布を外部から与える
C           2001/01/16  小高正嗣; ・散逸加熱(DTPVIS)を考慮
C                                 ・運動エネルギー収支を計算(EMONIT)
C**********************************************************************
C
      IMPLICIT REAL*8 ( A-H, O-Z )
C
C**********************************************************************
C
#include      "grid_size_M.f"
C
      PARAMETER ( NZPROF = 100 )!
      parameter ( NBND = 16 )   ! IR band number (15μm)
      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
C      parameter ( SIGMA  = 5.67D-8  )
      parameter ( PI = 3.14159263D0 )
C      parameter ( TSOL   = 5760.0D0 ) ! surface temp. of Sun
      parameter ( DAY    = 86400.0D0)
C      parameter ( HOUR   = 3600.0D0 ) 

C======================================================================
C
C *** Basic State Variables. ***
C
      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 )
C
      DIMENSION FCOLI( -2:NXB ),
     \          TSFC ( -2:NXB ),
     \          QVSFC( -2:NXB )
C
      DIMENSION
     \   FDZ  ( -2:NZB ), FDZM ( -2:NZB )
C
C *** Array for Pressure equation solver. ***
C
      DIMENSION
     \   VMODE ( -2:NZB, NZ  ), TMODE ( -2:NZB, NZ  ), VEIGEN (NZ),
     \   FMODE ( -2:NXB, NZ )
C
      DIMENSION
     \   TRIGS( NX+15 ), IFAC(15)
C
C
C *** Primary Time-dependent Variables. ***
C
      DIMENSION
     \     U ( -2:NXB, -2:NZB, NROT ),
     \     V ( -2:NXB, -2:NZB, NROT ),
     \     W ( -2:NXB, -2:NZB, NROT ),
     \  PPAI ( -2:NXB, -2:NZB, NROT ),
     \  TPOT ( -2:NXB, -2:NZB, NROT ),
     \  QVAP ( -2:NXB, -2:NZB, NROT ),
     \  QCLW ( -2:NXB, -2:NZB, NROT ),
     \  QRAI ( -2:NXB, -2:NZB, NROT ),
     \ CDTURB( -2:NXB, -2:NZB, NROT )

      DIMENSION
     \     CDRAGV( -2:NXB ),
     \     RICHD ( -2:NXB )
C
C *** Monitor variables. ***
C
      CHARACTER*10 NMZPRF(NZPROF)
      DIMENSION
     \   ZPROF ( -2:NZB, NZPROF )
      DIMENSION
     \   AXU    ( -2:NZB ), AXV    ( -2:NZB ), AXW    ( -2:NZB ),
     \   AXTPOT ( -2:NZB ), AXMSE  ( -2:NZB ),
     \   AXQVAP ( -2:NZB ), AXQCLW ( -2:NZB ), AXQRAI ( -2:NZB )
      DIMENSION
     \   IPQCLW ( -2:NXB ), IPQRAI ( -2:NXB )
      DIMENSION
     \   FUMOM ( -2:NXB,2 ), FVMOM ( -2:NXB,2 ),
     \   FHEAT ( -2:NXB,2 ), FEVAP ( -2:NXB,2 ), FWPREC( -2:NXB,2 )
      DIMENSION
     \   FUMOMA( -2:NXB,2 ), FVMOMA( -2:NXB,2 ),
     \   FHEATA( -2:NXB,2 ), FEVAPA( -2:NXB,2 ), FWPREA( -2:NXB,2 )
C
C
C *** for Cloudphysical Calculation. ***
C
      DIMENSION
     \     TPSTE ( 75:175, 0:NZ )

C
C *** for Radiative Transfer Calculation ***
C
      DIMENSION
     \   OPL   (-2:NZB)

      DIMENSION
     \   TAUN    (0:NZ+1,0:NZ+1,NBND)
c     \   TAUN_SR1(0:NZ+1,0:NZ+1,NBNDSR1),
c     \   TAUN_SR2(0:NZ+1,0:NZ+1,NBNDSR2),
c     \   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
     \   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
     \   FXRADD ( -2:NXB, -2:NZB ), SXRADD ( -2:NXB, -2:NZB ),
     \   FXQRADD( -2:NXB, -2:NZB ), SXQRADD( -2:NXB, -2:NZB ),
     \   SXSOL  ( -2:NXB, -2:NZB )

C
C *** for Surface Temprature Calculation ***
C
      dimension
     \   TEMPG(-2:NXB, NZG)

      dimension
     \   GGRID(-1:NZGB), GDZ(0:NZGB), GDZM(0:NZG)


C
C *** Physical Constants. ***
C
      NAMELIST /CONST1/ GRAV, PSD, AMDRY, RDRY, CPDRY, RMVAP, ALATNT
C
C *** Model Specification ***
C
      NAMELIST /MODEL/ NNX, NNZ, DDX
C
C *** Parameters to be specified in 'SETCS2'. ***
C
C      NAMELIST /CONST2/ CDRAG, VSFC0, PSFC, TSFC, QVSFC, FCOLI
      NAMELIST /CONST2/ CDRAG, VSFC0, PSFC
C
C *** Parameters to be specified in 'SETCS2'. ***
C
C=====================================================================
C
      DATA  ( NMZPRF(IST), IST = 11, 16 )
     \   / 'UMEAN', 'DUADV', 'DUVIS', 'DUNLV', 'DUCORI', 'DUFRIC' /
      DATA  ( NMZPRF(IST), IST = 21, 26 )
     \   / 'VMEAN', 'DVADV', 'DVVIS', 'DVNLV', 'DVCORI', 'DVFRIC' /
      DATA  ( NMZPRF(IST), IST = 31, 33 )
     \   / 'EKVIS', 'EKNLV', 'EKBUOY' /
      DATA  ( NMZPRF(IST), IST = 41, 50 )
     \   / 'TPMEAN', 'DTPADV', 'DTPDIF', 'DTPDNL', 'DTPAD0', 'DTPDI0',
     \     'DTPVIS', 'DTIRAD', 'DTSRAD', 'FHEAT' /
      DATA  ( NMZPRF(IST), IST = 51, 59 )
     \   / 'QVMEAN', 'DQVADV', 'DQVDIF', 'DQVDNL', 'DQVAD0', 'DQVDI0',
     \     'DQVNFL', 'DTQFAL', 'FEVAP' /
      DATA  ( NMZPRF(IST), IST = 61, 65 )
     \   / 'QCMEAN', 'DQCADV', 'DQCDIF', 'DQCDNL', 'DQCNFL' /
      DATA  ( NMZPRF(IST), IST = 71, 75 )
     \   / 'QRMEAN', 'DQRADV', 'DQRDIF', 'DQRDNL', 'DQRNFL' /
      DATA  ( NMZPRF(IST), IST = 81, 85 )
     \   / 'KSMEAN', 'DKSADV', 'DKSDIF', 'DKSDNL', 'DKSGEN' /
      DATA  ( NMZPRF(IST), IST = 91, 96 )
     \   / 'DQRFAL', 'PCCON',  'PVEVC',  'PVEVR',  'PCONV',  'PCOAG' /
C
C======================================================================
C
C     DO 9091 IRUN = 1, 10
C
c      write(0,*) 'IRUN ?'
C
C
      READ(5,*) IRUN
C
C
      WRITE(6,9601)
      WRITE(6,*) ' *** SERIES NUMBER', IRUN,' BEGINS HERE ***'
 9601 FORMAT ( 1H1 )
C
      CALL FOPEN ( IRUN )
C
C *** SET UP ***
C
      CALL SETCST
     O   ( GRAV, PSD, AMDRY, RDRY, CPDRY, RMVAP, ALATNT )
C
      CALL SETGRD
     O   ( ZGRID(-2), FDZ(-2), FDZM(-2), DX,
     D     NZ, NZB                           )
C
      CALL SETCS2
     O   ( CDRAG, VSFC0, TSFC(-2), QVSFC(-2),
     O     FCOLI,
     O     PSFC,
     C     RMVAP,
     D     NX, NXB                      )
C
C *** Dust Optical Paremeter ***(2000/11/09)
C 
      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 )
C
C *** For Surface Temp. Calculation Set Up ***(2000/10/31)
C 
      call SETCSG 
     I   ( DAY   ,
     O     ALBEDO, DIFG, CPGRND, DENSG, CDTURBG, SKDPS )

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

      DATA ( GGRID(IZ), IZ=-1,NZGB )
     1     / -6.0D0 ,
     2       -6.0D0 , -4.0D0 , -2.7D0 , -1.8D0, -1.2D0, 
     3       -0.79D0, -0.53D0, -0.35D0, -0.2D0, -0.1D0, 
     4        0.0D0 ,  0.0D0 /

      call SETGGD
     I   ( SKDPS, GGRID, 
     O     GDZ  , GDZM, 
     D     NZG  , NZGB       ) 
C
C
      WRITE(16, MODEL  )
      WRITE(16, CONST1 )
      WRITE(16, CONST2 )
C
C      write(0,*) 'i came here 1 '
C
      CALL VBASIC
     I   ( ZGRID, PSFC,
     O     PRES0(-2), DENS0(-2), TEMP0(-2), TPOT0(-2), PPAI0(-2),
     O     QVAP0(-2), RVAP0(-2),
     C     GRAV , RDRY , CPDRY, PSD,  RMVAP,
     D     NZ, NZB                              )
C
C      write(0,*) 'i came here 2 '
C
!-----------------------------------------------------------------      
!*** For Radiative Transfer Calculation Set Up***(2000/10/31)
! (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   / 

!-----------------------------------------------------------------      
!*** For CO2 15 μm band Set Up *** 

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

     
      DO 101 IZ = -2, NZB
         WRITE(6,601) IZ, ZGRID(IZ), PRES0(IZ), DENS0(IZ), TEMP0(IZ),
     \                    TPOT0(IZ), PPAI0(IZ), QVAP0(IZ), RVAP0(IZ)
         WRITE(16,*) IZ, ZGRID(IZ), PRES0(IZ), DENS0(IZ), TEMP0(IZ),
     \                     TPOT0(IZ), PPAI0(IZ), QVAP0(IZ), RVAP0(IZ)
  101 CONTINUE
C
  601 FORMAT ( I4, 8E15.7)
C
      call PSETUP
     B   (  DENS0,
     G      FDZ  , FDZM,
     M      VMODE, TMODE, VEIGEN, FMODE, TRIGS, IFAC )
C
C--------------------------------------------------------------------
C
      NNX = NX
      NNZ = NZ
      DDX = DX
C
      WRITE(6, MODEL  )
      WRITE(6, CONST1 )
C
C====================================================================
C *** SPECIFICATION OF TIME INTEGRATION ***
C TIME0    start model time (sec)
C NLOOP1   the number of time step every output
C NLOOP2   the number of data output 
C DTIME    time interval (sec)
C TIMES    offset time from LT=00:00 (sec)
C
      WRITE(6,*)'TIME0, NLOOP1, NLOOP2, DTIME, NTMONI, NTKUBU, TIMES'
C
C
      READ (5,*) TIME0, NLOOP1, NLOOP2, DTIME, NTMONI , NTKUBU, TIMES
C
C
      WRITE(6,*) TIME0, NLOOP1, NLOOP2, DTIME, NTMONI , NTKUBU, TIMES
      WRITE(16,*) TIME0, NLOOP1, NLOOP2, DTIME, NTMONI , NTKUBU, TIMES
C
C YLAT     latitude parameter
C ALS      Mars Ls parameter
C
      WRITE(6,*)'YLAT, ALS'
      READ (5,*) YLAT, ALS
      WRITE(6,*) YLAT, ALS

      YLAT = YLAT * PI/180.0D0  ! transfer from degree to rad
      ALS  = ALS  * PI/180.0D0  ! transfer from degree to rad
C
C--------------------------------------------------------------------
C *** INITIAL CONDITIONS ***
C--------------------------------------------------------------------

C
      IF ( TIME0  .EQ.  0.0D0 ) THEN
C
      CALL DAINIT
     I   (  1    ,
     O      TIME,
     O      U, V, W, TPOT, QVAP, QCLW, QRAI, CDTURB,
     B      PRES0, DENS0, TEMP0, TPOT0, PPAI0, QVAP0, RVAP0, ZGRID,
     C      RMVAP,
     D      NX, NZ, NXB, NZB, NROT           )
C
C      do IZ = 1, NZG
C         do IX = -2, NZB
C            TEMPG(IX,IZ) = 220.0D0
C         end do
C      end do

      OPEN(11,FILE='1DGRAND.DATA')
      READ(11,*) ( TEMPG(-2, IZ), IZ=1,NZG )
      CLOSE(11)

      do IZ = 1, NZG
         do IX = -1, NXB
            TEMPG(IX,IZ) = TEMPG(-2, IZ)
         end do
      end do

      ELSE
C
      CALL FLREAD
     I   (  TIME0,  1,
     O      TIME,
     O      U, V, W, TPOT, QVAP, QCLW, QRAI, PPAI, CDTURB, TEMPG,
     D      NX, NZ, NXB, NZB, NZG, NROT           )
C
      END IF

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

c      write(0,*) ' came here 0001 '

C
      IF ( IRUN .EQ. 1 ) then
         CALL FLWRIT
     I      (  TIME ,  1,
     I         U, V, W, TPOT, QVAP, QCLW, QRAI, PPAI, CDTURB, TEMPG,
     D         NX, NZ, NXB, NZB, NZG, NROT           )
      END IF
C
C----------------------------------------------------------------------
C *** TIME INTEGRATION ***
C----------------------------------------------------------------------
C
      NTLAG = 20
C
C
      DO 9999 INL1 = 1, NLOOP1
          WRITE(6,*) INL1, TIME
      DO 999 INL2 = 1, NLOOP2
C         WRITE(0,*) INL1, INL2, TIME
C
      ITL = INL2 + NLOOP2*(INL1-1)
C
C     デバッグ用出力を NDP 回毎に行う。
C
c     NDP = 1
c     IDUMP = MOD( (ITL-1), NDP )
c     IF ( IDUMP .EQ. 0 ) THEN
c        IFDUMP = 1
c     ELSE
         IFDUMP = 0
c     END IF
C
      ILAG = MOD( (ITL-1), NTLAG )
      IF ( ILAG .EQ. 0 ) THEN
C
      IT1 = 1
      IT2 = 1
      IC1 = 1
      IC2 = 2
      DTI2 = DTIME / 2.0D0
C
C
      CALL TSTEP
     I   ( IT1, IT2, IC1, IC2, DTI2,
     I     TIME, IFDUMP, TIMES, TIME0, 
     C     DX, GRAV, CPDRY, ALATNT, RMVAP, 
     C     FCOLI, CDRAG, VSFC0, TSFC, QVSFC,
     C     PSFC, RDRY,
     V     U   , V   , W   , PPAI, TPOT, QVAP, QCLW, QRAI, CDTURB,
     V     CDRAGV, RICHD , TEMPG , 
     S     DTPRAD, DTIRAD, DTSRAD, DTQIRAD, DTQSRAD,
     S     FXRADD, SXRADD, FXQRADD,SXQRADD, SXSOL  ,
     F     FUMOM , FVMOM , FHEAT , FEVAP , FWPREC,
     F     FUMOMA, FVMOMA, FHEATA, FEVAPA, FWPREA,
     B     ZGRID, FDZ  , FDZM , TEMP0, TPOT0,
     B     DENS0, PPAI0, QVAP0, pres0,
     B     GDZ  , GDZM ,
     B     TAUN , OPL  , 
     C     WNU  ,
     C     WNU_SR1, SNU_SR1, RNU_SR1,
     C     WNU_SR2, SNU_SR2, RNU_SR2,
     C     WNU_SR3, SNU_SR3, RNU_SR3,
     C     QEXT_SR, QEXT_IR1, QEXT_IR2, 
     C     QALB_SR, QALB_IR1, QALB_IR2, 
     C     QASM_SR, QASM_IR1, QASM_IR2,
     C     QALB_SRE, QASM_SRE, QDENS  , QREFF   , QRMOD  ,  
     C     S00   , ALS   , YLAT   , EXCTL , OBLIC , OMGA0, 
     C     ALBEDO, CPGRND, CDTURBG, DENSG ,
     M     VMODE, TMODE, FMODE, TRIGS, IFAC,
     P     TPSTE ,
     B     CTHKIN, CLATHR, CPWKIN,
     B     SETHER, SELATN, SEPOTW,
     B     DEKVIS, DEKNLV, DPWFAL,
     B     WCVTOL, WSVAPR, WSLIQD, XQFILL,
     B     ZPROF                               )
C
C
C
      ELSE
C
      IT1 = MOD( (ITL-1), 2 ) + 1
      IT2 = MOD(  ITL   , 2 ) + 1
      IC1 = IT2
      IC2 = IT2
      DTI = DTIME
C
C
      CALL TSTEP
     I   ( IT1, IT2, IC1, IC2, DTI,
     I     TIME, IFDUMP, TIMES, TIME0, 
     C     DX, GRAV, CPDRY, ALATNT, RMVAP,
     C     FCOLI, CDRAG, VSFC0, TSFC, QVSFC,
     C     PSFC, RDRY,
     V     U   , V   , W   , PPAI, TPOT, QVAP, QCLW, QRAI, CDTURB,
     V     CDRAGV, RICHD , TEMPG , 
     S     DTPRAD, DTIRAD, DTSRAD, DTQIRAD, DTQSRAD,
     S     FXRADD, SXRADD, FXQRADD,SXQRADD, SXSOL  ,
     F     FUMOM , FVMOM , FHEAT , FEVAP , FWPREC,
     F     FUMOMA, FVMOMA, FHEATA, FEVAPA, FWPREA,
     B     ZGRID, FDZ  , FDZM , TEMP0, TPOT0,
     B     DENS0, PPAI0, QVAP0, pres0,
     B     GDZ  , GDZM ,
     B     TAUN,  OPL  , 
     C     WNU  ,
     C     WNU_SR1, SNU_SR1, RNU_SR1,
     C     WNU_SR2, SNU_SR2, RNU_SR2,
     C     WNU_SR3, SNU_SR3, RNU_SR3,
     C     QEXT_SR, QEXT_IR1, QEXT_IR2, 
     C     QALB_SR, QALB_IR1, QALB_IR2, 
     C     QASM_SR, QASM_IR1, QASM_IR2,
     C     QALB_SRE, QASM_SRE, QDENS  , QREFF   , QRMOD  ,  
     C     S00   , ALS   , YLAT   , EXCTL , OBLIC , OMGA0, 
     C     ALBEDO, CPGRND, CDTURBG, DENSG ,
     M     VMODE, TMODE, FMODE, TRIGS, IFAC,
     P     TPSTE,
     B     CTHKIN, CLATHR, CPWKIN,
     B     SETHER, SELATN, SEPOTW,
     B     DEKVIS, DEKNLV, DPWFAL,
     B     WCVTOL, WSVAPR, WSLIQD, XQFILL,
     B     ZPROF                               )
C
C
      END IF
C
      TIME = TIME + DTIME / 2.0D0
C
      IKB = MOD ( ITL, NTKUBU )
      IF ( IKB .EQ. 0 ) THEN
         CALL KUBUN
     I      (  QCLW, QRAI, IC2,
     O         IPQCLW, IPQRAI,
     B         DENS0,
     D         NX, NZ, NXB, NZB, NROT )
         WRITE(71) TIME, IPQCLW, IPQRAI
      END IF
C
      IMN = MOD ( ITL, NTMONI )
      IF ( IMN .EQ. 0 ) THEN
C
C      WRITE(61) TIME, NMZPRF, ZPROF
      WRITE(61) TIME, ZPROF
      WRITE(72) TIME, FUMOMA, FVMOMA, FHEATA, FEVAPA, FWPREA
C
      DO 201 IST = 1, NZPROF
      DO 201 IZ = -2, NZB
         ZPROF (IZ,IST) = 0.0D0
  201 CONTINUE
C
      DO 202 IX = 1, NX
         FUMOMA(IX,1) = 0.0D0
         FVMOMA(IX,1) = 0.0D0
         FHEATA(IX,1) = 0.0D0
         FEVAPA(IX,1) = 0.0D0
         FWPREA(IX,1) = 0.0D0
         FUMOMA(IX,2) = 0.0D0
         FVMOMA(IX,2) = 0.0D0
         FHEATA(IX,2) = 0.0D0
         FEVAPA(IX,2) = 0.0D0
         FWPREA(IX,2) = 0.0D0
  202 CONTINUE
C
      CALL MONIT
     I   ( U     , V     , W     , PPAI ,
     I     TPOT  , QVAP  , QCLW  , QRAI ,
     I     IC2   , TIME  ,
     O     AXU   , AXV   , AXW   ,
     O     AXTPOT, AXQVAP, AXQCLW, AXQRAI, AXMSE ,
     O     EKINET, ETHERM, ELATNT, EPOTW ,
     M     CTHKIN, CLATHR, CPWKIN,
     M     SETHER, SELATN, SEPOTW,
     M     DEKVIS, DEKNLV, DPWFAL,
     B     DENS0 , TPOT0 , TEMP0 , PPAI0 , QVAP0 , ZGRID,
     C     GRAV  , CPDRY , RMVAP , ALATNT,
     G     DX    , FDZ   , FDZM  ,
     D     NX    , NZ    , NXB   , NZB   ,  NROT   )
C
      CALL MONITW
     O   ( WATOTL, WAVAPR, WALIQD,
     M     WCVTOL, WSVAPR, WSLIQD, XQFILL,
     I     QVAP  , QCLW  , QRAI  ,
     I     IC2   , TIME  ,
     B     DENS0 , FDZ   , QVAP0 ,
     G     DX    ,
     D     NX, NZ, NXB, NZB, NROT )
C
      END IF
C
  999 CONTINUE
C
      IWRIT = IC2
C
      CALL FLWRIT
     I   (  TIME ,  IWRIT,
     I      U, V, W, TPOT, QVAP, QCLW, QRAI, PPAI, CDTURB, TEMPG,
     D      NX, NZ, NXB, NZB, NZG, NROT           )
C
C
 9999 CONTINUE
C
C
      CALL FCLOSE
C
C 9091 CONTINUE
C
      STOP
      END
      SUBROUTINE TSTEP
     I   ( IT1, IT2, IC1, IC2, DTI,
     I     TIME, IFDUMP, TIMES, TIME0, 
     C     DX, GRAV, CPDRY, ALATNT, RMVAP,
     C     FCOLI, CDRAG, VSFC0, TSFC, QVSFC,
     C     PSFC, RDRY,
     V     U   , V   , W   , PPAI, TPOT, QVAP, QCLW, QRAI, CDTURB,
     V     CDRAGV, RICHD, TEMPG, 
     S     DTPRAD, DTIRAD, DTSRAD, DTQIRAD, DTQSRAD,
     S     FXRADD, SXRADD, FXQRADD,SXQRADD, SXSOL  ,
     F     FUMOM , FVMOM , FHEAT , FEVAP , FWPREC,
     F     FUMOMA, FVMOMA, FHEATA, FEVAPA, FWPREA,
     B     ZGRID, FDZ  , FDZM , TEMP0, TPOT0,
     B     DENS0, PPAI0, QVAP0, pres0,
     B     GDZ  , GDZM ,
     B     TAUN,  OPL  ,
     C     WNU  ,
     C     WNU_SR1, SNU_SR1, RNU_SR1,
     C     WNU_SR2, SNU_SR2, RNU_SR2,
     C     WNU_SR3, SNU_SR3, RNU_SR3,
     C     QEXT_SR, QEXT_IR1, QEXT_IR2, 
     C     QALB_SR, QALB_IR1, QALB_IR2, 
     C     QASM_SR, QASM_IR1, QASM_IR2,
     C     QALB_SRE, QASM_SRE, QDENS  , QREFF   , QRMOD  ,  
     C     S00   , ALS   , YLAT   , EXCTL , OBLIC , OMGA0, 
     C     ALBEDO, CPGRND, CDTURBG, DENSG ,
     M     VMODE, TMODE, FMODE, TRIGS, IFAC,
     P     TPSTE ,
     B     CTHKIN, CLATHR, CPWKIN,
     B     SETHER, SELATN, SEPOTW,
     B     DEKVIS, DEKNLV, DPWFAL,
     B     WCVTOL, WSVAPR, WSLIQD, XQFILL,
     B     ZPROF                               )
C
C*********************************************************************
      IMPLICIT REAL*8 ( A-H, O-Z )
C*********************************************************************
C
#include      "grid_size_M.f"
C
      PARAMETER ( NZPROF = 100 )
      parameter ( NBND = 16 )   ! IR band number (15μm)
      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)
C      parameter ( HOUR   = 3600.0D0 ) 

C
C======================================================================
C
C *** Basic State Variables. ***
C
      DIMENSION
     \   ZGRID ( -2:NZB ),
     \   TEMP0 ( -2:NZB ),
     \   TPOT0 ( -2:NZB ),
     \   DENS0 ( -2:NZB ),
     \   PPAI0 ( -2:NZB ),
     \   Pres0 ( -2:NZB ),
     \   QVAP0 ( -2:NZB ),
     \   QADVM ( -2:NZB )
C
      DIMENSION  CD00 ( -2:NZB )
C
      DIMENSION
     \   FCOLI ( -2:NXB ),
     \   TSFC  ( -2:NXB ),
     \   QVSFC ( -2:NXB )
C
      DIMENSION
     \   FDZ  ( -2:NZB ), FDZM ( -2:NZB )
C
C *** Array for Pressure equation solver. ***
C
      DIMENSION
     \   VMODE ( -2:NZB, NZ  ), TMODE ( -2:NZB, NZ  ),
     \   FMODE ( -2:NXB, NZ )

c     dimension  UUFT(-2:NXB,NZ), ANFT(-2:NXB,NZ)
      dimension  UUFT(NZ,NX), ANFT(NZ,NX)
c
      dimension
     \   TRIGS ( NX+15 ), IFAC( 15 )
C
C *** Primary Time-dependent Variables. ***
C
      DIMENSION
     \     U ( -2:NXB, -2:NZB, NROT ),
     \     V ( -2:NXB, -2:NZB, NROT ),
     \     W ( -2:NXB, -2:NZB, NROT ),
     \  PPAI ( -2:NXB, -2:NZB, NROT ),
     \  TPOT ( -2:NXB, -2:NZB, NROT ),
     \  QVAP ( -2:NXB, -2:NZB, NROT ),
     \  QCLW ( -2:NXB, -2:NZB, NROT ),
     \  QRAI ( -2:NXB, -2:NZB, NROT ),
     \ CDTURB( -2:NXB, -2:NZB, NROT )

      DIMENSION
     \     CDRAGV( -2:NXB ),
     \     RICHD ( -2:NXB )

C *** Surface Temprature ***

      dimension
     \   TEMPG(-2:NXB, NZG)

      dimension
     \   GDZ(0:NZGB), GDZM(0:NZG)

C
C *** Monitor variables ***
C
      DIMENSION
     \   ZPROF( -2:NZB, NZPROF )
C
C *** Secondary Variables. ***
C
      DIMENSION
     \     UFLX ( -2:NXB, -2:NZB ),
     \     WFLX ( -2:NXB, -2:NZB ),
     \     UM1  ( -2:NXB, -2:NZB ),
     \     VM1  ( -2:NXB, -2:NZB ),
     \     WM1  ( -2:NXB, -2:NZB ),
     \     UM2  ( -2:NXB, -2:NZB ),
     \     VM2  ( -2:NXB, -2:NZB ),
     \     WM2  ( -2:NXB, -2:NZB )
C
      DIMENSION
     \     PFUNC ( -2:NXB, -2:NZB ),
     \     BUOY  ( -2:NXB, -2:NZB ),
     \     BUOYW ( -2:NXB, -2:NZB )
C
      DIMENSION
     \     CONVM ( -2:NXB, -2:NZB ),
     \     EDFM1 ( -2:NXB, -2:NZB ),
     \     EDFM2 ( -2:NXB, -2:NZB ),
     \     STABZ ( -2:NXB, -2:NZB ),
     \     CD1   ( -2:NXB, -2:NZB ),
     \     CD2   ( -2:NXB, -2:NZB )
C
C *** Surface Fluxes. ***
C
      DIMENSION
     \   FUMOM ( -2:NXB,2 ), FVMOM ( -2:NXB,2 ),
     \   FHEAT ( -2:NXB,2 ), FEVAP ( -2:NXB,2 ), FWPREC( -2:NXB,2 )
C
      DIMENSION
     \   FUMOMA( -2:NXB,2 ), FVMOMA( -2:NXB,2 ),
     \   FHEATA( -2:NXB,2 ), FEVAPA( -2:NXB,2 ), FWPREA( -2:NXB,2 )
C
C *** Radiative Cooling. ***
C
      DIMENSION
     \   DTPRAD ( -2:NXB, -2:NZB ),
     \   TPMEAN ( -2:NZB ), ALMEAN( -2:NZB ),
     \   ALDEVI ( -2:NZB ), COOLZ ( -2:NZB )

      DIMENSION
     \   DTIRAD ( -2:NXB, -2:NZB ),
     \   DTSRAD ( -2:NXB, -2:NZB ),
     \   DTQIRAD( -2:NXB, -2:NZB ),
     \   DTQSRAD( -2:NXB, -2:NZB )
C
C *** Disspative Heating. ***
C
      DIMENSION
     \     DTPVIS ( -2:NXB, -2:NZB )
C
C *** Time Derivatives. ***
C
      DIMENSION
     \     DUTOTL( -2:NXB, -2:NZB ), DVTOTL( -2:NXB, -2:NZB ),
     \     DWTOTL( -2:NXB, -2:NZB )
      DIMENSION
     \     DUADV ( -2:NXB, -2:NZB ), DUVIS ( -2:NXB, -2:NZB ),
     \     DUNLV ( -2:NXB, -2:NZB ),
     \     DVADV ( -2:NXB, -2:NZB ), DVVIS ( -2:NXB, -2:NZB ),
     \     DVNLV ( -2:NXB, -2:NZB ),
     \     DWADV ( -2:NXB, -2:NZB ), DWVIS ( -2:NXB, -2:NZB ),
     \     DWNLV ( -2:NXB, -2:NZB ),
     \     DUCORI( -2:NXB, -2:NZB ), DVCORI( -2:NXB, -2:NZB )
C
      DIMENSION
     \     DUFRIC ( -2:NXB,2 )     , DVFRIC( -2:NXB,2 )
C
      DIMENSION
     \     DTQADV ( -2:NXB, -2:NZB ), DTQDIF ( -2:NXB, -2:NZB ),
     \     DTQAD0 ( -2:NXB, -2:NZB ), DTQDI0 ( -2:NXB, -2:NZB ),
     \     DTQDNL ( -2:NXB, -2:NZB ), DTQFAL ( -2:NXB, -2:NZB )
C
      DIMENSION
     \     DQFILL ( -2:NXB, -2:NZB ), QSUMPN ( -2:NXB, -2:NZB )
C
C *** for the kinetic energy buget. (2001/05/01) ***
C
      dimension
     \     EKVIS ( -2:NXB, -2:NZB ),
     \     EKNLV ( -2:NXB, -2:NZB ),
     \     EKBUOY( -2:NXB, -2:NZB )
C
C *** for the calculation of turbulent diffusivity. ***
C
      DIMENSION
     \     SSTURB ( -2:NXB, -2:NZB )
C
C *** for Cloudphysical Calculation. ***
C
      DIMENSION
     \     TPSTE ( -25:75, 0:NZ )
C
      DIMENSION
     \     TPOT1 ( -2:NXB, -2:NZB ),
     \     QVAP1 ( -2:NXB, -2:NZB ),
     \     QCLW1 ( -2:NXB, -2:NZB ),
     \     QRAI1 ( -2:NXB, -2:NZB )
C
      DIMENSION
     \     PRMIC ( -2:NXB, -2:NZB ),
     \     PCONV ( -2:NXB, -2:NZB ),
     \     PCOAG ( -2:NXB, -2:NZB ),
     \     PCCON ( -2:NXB, -2:NZB ),
     \     PVEVC ( -2:NXB, -2:NZB ),
     \     PVEVR ( -2:NXB, -2:NZB )
C
C
      DIMENSION
     \     VRFALL( -2:NXB, -2:NZB ),
     \     FXRAIN( -2:NXB, -2:NZB )
C
C *** for Radiative Transfer Calculation ***
C
      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),
     \   WNU(NBND),
     \   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
     \   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
     \   TEMP  ( -2:NXB, -2:NZB ), BBRAD ( -2:NXB, -2:NZB ),
     \   BBRAD_S( -2:NXB )       , FXRADU( -2:NXB, -2:NZB ),
     \   FXRADD( -2:NXB, -2:NZB ), FXRADN( -2:NXB, -2:NZB )

      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) 

C
C *** For Scattering Matrix Array *** (2000/11/09)
C
      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)
C
C *** Work space. ***
C
      DIMENSION
     \     FLQX ( -2:NXB, -2:NZB ),
     \     FLQZ ( -2:NXB, -2:NZB )
      DIMENSION
     \     QLAP  ( -2:NXB, -2:NZB ),
     \     CDQNLX( -2:NXB, -2:NZB ),
     \     CDQNLZ( -2:NXB, -2:NZB ),
     \     DIVMU ( -2:NXB, -2:NZB ),
     \     DIVMW ( -2:NXB, -2:NZB ),
     \     DIV   ( -2:NXB, -2:NZB ),
     \     PDIV  ( -2:NXB, -2:NZB ),
     \     PDIVM ( -2:NXB, -2:NZB ),
     \     PFUNCM( -2:NXB, -2:NZB ),
     \     ALPHA ( -2:NXB, -2:NZB ),
     \     GAMMA ( -2:NXB, -2:NZB )
C
C=====================================================================
C
C *** WORK AREA SPECIFICATION ***
C
C----------------------------------------------------------------------
C
      COMMON /WORK/ WORK( -2:NXB, -2:NZB, 17 )
C
C
      EQUIVALENCE
     \ ( WORK(-2,-2, 1), UM1  ,  DUVIS,           UFLX, PRMIC     ),
     \ ( WORK(-2,-2, 2), VM1  ,  DVVIS,           WFLX, PCONV     ),
     \ ( WORK(-2,-2, 3), WM1  ,  DWVIS, DQFILL,   FLQX, PCOAG     ),
     \ ( WORK(-2,-2, 4), UM2  ,  DUNLV, QSUMPN,   FLQZ, PCCON     ),
     \ ( WORK(-2,-2, 5), VM2  ,  DVNLV,         DTQADV, PVEVC     ),
     \ ( WORK(-2,-2, 6), WM2  ,  DWNLV,         DTQAD0, PVEVR     ),
     \ ( WORK(-2,-2, 7), DUADV, PFUNCM,         DTQDIF            ),
     \ ( WORK(-2,-2, 8), DVADV,  PDIVM,         DTQDI0, VRFALL    ),
     \ ( WORK(-2,-2, 9), DWADV,  PDIV ,         DTQDNL, FXRAIN    ),
     \ ( WORK(-2,-2,10), DIVMU,  CD1                              ),
     \ ( WORK(-2,-2,11), DIVMW,  CD2  , DUCORI,  CONVM,  TPOT1    ),
     \ ( WORK(-2,-2,12), EDFM1,         DVCORI,  ALPHA,  QVAP1    ),
     \ ( WORK(-2,-2,13), EDFM2,          BUOY ,  GAMMA,  QCLW1    ),
     \ ( WORK(-2,-2,14), STABZ,  DIV  ,  BUOYW,          QRAI1    )
C     \ ( WORK(-2,-2,14), STABZ,  DIV  ,  BUOYW, DTPRAD,  QRAI1    )
C    \ ( WORK(-2,-2,15), QLAP                                     ),
C    \ ( WORK(-2,-2,16), CDQNLX                                   ),
C    \ ( WORK(-2,-2,17), CDQNLZ                                   )
C
C====================================================================
C


      IF ( IT1 .EQ. IT2 ) THEN
         DTBGT = DTI
      ELSE
         DTBGT = DTI / 2.0D0
      END IF
C
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
      CALL CLWMID
     I   ( U, V, W,
     I     IT1,
     O     UM1, UM2, VM1, VM2, WM1, WM2,
     D     NX, NZ, NXB, NZB, NROT )
C
      CALL CLWADV
     I   ( UM1, VM1, WM1,
     I     UM2, VM2, WM2,
     O     DUADV, DVADV, DWADV,
     W     DIVMU, DIVMW,
     B     DENS0, FDZ, FDZM,
     G     DX,
     D     NX, NZ, NXB, NZB      )
C
C
      CALL CDCOEF
     I   ( CDTURB,
     I     IT1   ,
     O     CD1   , CD2   ,
     G     ZGRID ,
     D     NX    , NZ    , NXB   , NZB   , NROT  )
C
      CALL ADUMP ( CD1 , 'CD1 ', NXB, NZB, 1, 1 , IFDUMP,TIME )
C     CALL ADUMP ( CD2 , 'CD2 ', NXB, NZB, 1, 1 , IFDUMP,TIME )
C
c --- Time Dependent TSFC----------------------------------------

C      DO IX = -2, NX+2
c         TSFC(IX) = 180.0D0 + 1.0D+2 * 
c     \   MAX( sin( 2.0d0 * PI * (TIME+2.88D+4) / 8.64D+4 ), 0.0D0 )
C         TSFC(IX) = 180.0D0 + 1.0D+2 * 
C     \              MAX( sin( 2.0d0 * PI * TIME / 8.64D+4 ), 0.0D0 )
C      END DO

c ---------------------------------------------------------------

      CALL CLDRAG
     I   ( U    , V    , TPOT, 
     I     IT2, 
     I     TPOT0, 
     O     RICHD, CDRAGV, 
     G     ZGRID, 
     C     GRAV , TSFC ,
     D     NX, NXB, NZB, NROT )

      CALL FXSURF3
     I   ( U,     V,     TPOT , QVAP ,
     I     IT2  , CDRAGV,
     O     FUMOM, FVMOM, FHEAT, FEVAP,
     B     DENS0, TPOT0, QVAP0, PPAI0,
     C     CPDRY,
     S     VSFC0,
     S     TSFC, QVSFC,
     D     NX, NZ, NXB, NZB, NROT    )
C
C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C     Heat and Water budget ( Surface flux terms )
C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
      DO 901 IX = 1, NX
         SETHER = SETHER + FHEAT(IX,1)          * DTBGT * DX
     \                   + FHEAT(IX,2)          * DTBGT * DX
         SELATN = SELATN + ALATNT * FEVAP(IX,1) * DTBGT * DX
     \                   + ALATNT * FEVAP(IX,2) * DTBGT * DX
         WSVAPR = WSVAPR + FEVAP(IX,1)          * DTBGT * DX
     \                   + FEVAP(IX,2)          * DTBGT * DX
         FHEATA(IX,1) = FHEATA(IX,1) + FHEAT(IX,1) * DTBGT
         FEVAPA(IX,1) = FEVAPA(IX,1) + FEVAP(IX,1) * DTBGT
         FUMOMA(IX,1) = FUMOMA(IX,1) + FUMOM(IX,1) * DTBGT
         FVMOMA(IX,1) = FVMOMA(IX,1) + FVMOM(IX,1) * DTBGT
         FHEATA(IX,2) = FHEATA(IX,2) + FHEAT(IX,2) * DTBGT
         FEVAPA(IX,2) = FEVAPA(IX,2) + FEVAP(IX,2) * DTBGT
         FUMOMA(IX,2) = FUMOMA(IX,2) + FUMOM(IX,2) * DTBGT
         FVMOMA(IX,2) = FVMOMA(IX,2) + FVMOM(IX,2) * DTBGT
  901 CONTINUE
C
      CALL CLVISC
     I   (   U,   V,   W,
     I     CD1, CD2, EDFM1,
     I     IT2,
     W     DIV,
     O     DUVIS, DVVIS, DWVIS,
     B     ZGRID, DENS0, FDZ, FDZM, CD00,
     G     DX,
     D     NX, NZ, NXB, NZB, NROT  )
C
C     CALL ADUMP ( DIV   , 'DIV  ', NXB, NZB, 1, 1 , IFDUMP,TIME )
      CALL ADUMP ( DUVIS , 'DUVIS', NXB, NZB, 1, 1 , IFDUMP,TIME )
C     CALL ADUMP ( DVVIS , 'DVVIS', NXB, NZB, 1, 1 , IFDUMP,TIME )
      CALL ADUMP ( DWVIS , 'DWVIS', NXB, NZB, 1, 1 , IFDUMP,TIME )
      CALL ADUMP ( DUADV , 'DUADV', NXB, NZB, 1, 1 , IFDUMP,TIME )
C     CALL ADUMP ( DVADV , 'DVADV', NXB, NZB, 1, 1 , IFDUMP,TIME )
      CALL ADUMP ( DWADV , 'DWADV', NXB, NZB, 1, 1 , IFDUMP,TIME )
C
C +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C     Kinetic energy dissipation due to viscosity.
C +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
      DO 2000 IZ = 0, NZ-1
         DMZ = DTBGT * FDZ(IZ) * DX**2 * DENS0(IZ)
      DO 2000 IX = 1, NX
         DEKVIS = DEKVIS +
     \      DMZ * (   U(IX,IZ,IT2) * DUVIS(IX,IZ)
     \              + V(IX,IZ,IT2) * DVVIS(IX,IZ) )
 2000 CONTINUE
C
      DO 2010 IZ = 1, NZ-1
         DMZ = DTBGT*FDZM(IZ)* DX**2 * (DENS0(IZ)+DENS0(IZ-1))/2.0D0
      DO 2010 IX = 1, NX
         DEKVIS = DEKVIS +
     \      DMZ * W(IX,IZ,IT2) * DWVIS(IX,IZ)
 2010 CONTINUE
C
C
      DO 2030 IX = 1, NX
         DEKVIS = DEKVIS +
     \    DTBGT * DX**2 *
     \    (  FDZ(0) * DENS0(0)
     \     * (   U(IX,0,IT2)*DUFRIC(IX,1) 
     \         + V(IX,0,IT2)*DVFRIC(IX,1) )
     \     + FDZ(NZ-1) * DENS0(NZ-1)
     \     * (   U(IX,NZ-1,IT2)*DUFRIC(IX,2)
     \         + V(IX,NZ-1,IT2)*DVFRIC(IX,2) ) )

 2030 CONTINUE

C
C +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
      SCVNLV = 3.0D0
C
      CALL NLVISC
     I   ( U, V, W, EDFM1, SCVNLV,
     I     IT2,
     O     DUNLV, DVNLV, DWNLV,
     B     DENS0, FDZ, FDZM,
     D     NX, NZ, NXB, NZB, NROT )
C
      CALL ADUMP ( DUNLV , 'DUNLV', NXB, NZB, 1, 1 , IFDUMP,TIME )
C     CALL ADUMP ( DVNLV , 'DVNLV', NXB, NZB, 1, 1 , IFDUMP,TIME )
      CALL ADUMP ( DWNLV , 'DWNLV', NXB, NZB, 1, 1 , IFDUMP,TIME )
C
C +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C     Kinetic energy dissipation due to computational viscosity.
C +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
      DO 3000 IZ = 0, NZ-1
         DMZ = DTBGT * FDZ(IZ) * DX**2 * DENS0(IZ)
      DO 3000 IX = 1, NX
         DEKNLV = DEKNLV +
     \      DMZ * (   U(IX,IZ,IT2) * DUNLV(IX,IZ)
     \              + V(IX,IZ,IT2) * DVNLV(IX,IZ) )
 3000 CONTINUE
C
      DO 3010 IZ = 1, NZ-1
         DMZ = DTBGT*FDZM(IZ)* DX**2 * (DENS0(IZ)+DENS0(IZ-1))/2.0D0
      DO 3010 IX = 1, NX
         DEKNLV = DEKNLV +
     \      DMZ * W(IX,IZ,IT2) * DWNLV(IX,IZ)
 3010 CONTINUE
C
C---------------------------------------------------------------------
C     Terms in the equations of motion.
C---------------------------------------------------------------------
C
      CALL CLBUOY
     I  (  TPOT(-2,-2,1), QVAP(-2,-2,1),
     I     QCLW(-2,-2,1), QRAI(-2,-2,1),
     I     IT1,
     O     BUOY(-2,-2), BUOYW(-2,-2),
     B     TPOT0(-2),
     C     GRAV, RMVAP,
     G     FDZ(-2),
     D     NX, NZ, NXB, NZB, NROT )
C
C+++++ Conversion to Kinetic energy ++++++++++++++++++++++++++++++++
C
      DO 1000 IZ = 1, NZ-1
         DMZT = FDZM(IZ) * DX**2 * (DENS0(IZ)+DENS0(IZ-1)) / 2.0D0
     \                   * DTBGT
         DO 1001 IX = 1, NX
            CTHKIN = CTHKIN + BUOY(IX,IZ) * W(IX,IZ,IT1)   * DMZT
 1001    CONTINUE
C
         DO 1002 IX = 1, NX
            CPWKIN = CPWKIN - BUOYW(IX,IZ) * W(IX,IZ,IT1)  * DMZT
 1002    CONTINUE
 1000 CONTINUE
C
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C     CALL ADUMP ( BUOY, 'BUOY', NXB, NZB, 1, 1 , IFDUMP,TIME )
C     WRITE(6,*) 'CLBUOY END'
C
      CALL CLCORI
     I   ( U, V,  IT1,
     O     DUCORI, DVCORI,
     I     FCOLI, ZGRID,
     D     NX, NZ, NXB, NZB, NROT )
C
C     CALL ADUMP ( DUCORI , 'DUCORI', NXB, NZB, 1, 1 , IFDUMP,TIME )
C     CALL ADUMP ( DVCORI , 'DVCORI', NXB, NZB, 1, 1 , IFDUMP,TIME )
C
C--------------------------------------------------------------------
C     Pressure diagnosis.
C--------------------------------------------------------------------
C
      DO 11 IX = -2, NX+2
         DUFRIC(IX,1) = FUMOM(IX,1) / DENS0(0) / DX / FDZ(0)
         DVFRIC(IX,1) = FVMOM(IX,1) / DENS0(0) / DX / FDZ(0)
         DUFRIC(IX,2) = FUMOM(IX,2) / DENS0(NZ-1) / DX / FDZ(NZ-1)
         DVFRIC(IX,2) = FVMOM(IX,2) / DENS0(NZ-1) / DX / FDZ(NZ-1)
   11 CONTINUE
C
C ---------------------------------------------------------------------
C  Time derivative of U, V and W except for the pressure gradient term.
C ---------------------------------------------------------------------
C
      DO 21 IZ = 0, NZ
      DO 21 IX = -2, NXB
         DUTOTL(IX,IZ) =   DUADV (IX,IZ) + DUVIS(IX,IZ)
     \                   + DUCORI(IX,IZ) + DUNLV(IX,IZ)
         DVTOTL(IX,IZ) =   DVADV (IX,IZ) + DVVIS(IX,IZ)
     \                   + DVCORI(IX,IZ) + DVNLV(IX,IZ)
         DWTOTL(IX,IZ) =   DWADV (IX,IZ) + DWVIS(IX,IZ)
     \                                   + DWNLV(IX,IZ)
     \                   + BUOY (IX,IZ)  + BUOYW(IX,IZ)
   21 CONTINUE
C
      DO 22 IX = 1, NX
         DUTOTL(IX,0) = DUTOTL(IX,0) + DUFRIC(IX,1)
         DVTOTL(IX,0) = DVTOTL(IX,0) + DVFRIC(IX,1)
         DUTOTL(IX,NZ-1) = DUTOTL(IX,NZ-1) + DUFRIC(IX,2)
         DVTOTL(IX,NZ-1) = DVTOTL(IX,NZ-1) + DVFRIC(IX,2)
   22 CONTINUE
C
      CALL BOUND ( DUTOTL, NX, NZ, NXB, NZB, 1, 1, 1, 1 )
      CALL BOUND ( DVTOTL, NX, NZ, NXB, NZB, 1, 1, 1, 1 )
C
C +++ BUDGET CALCULATION +++
C
      CALL CLZPRF
     I (     U, IC1, 11, 0, 0.0D0, ZPROF,
     D       NX, NZ, NXB, NZB, NROT, NZPROF )
      CALL CLZPRF
     I (     V, IC1, 21, 0, 0.0D0, ZPROF,
     D       NX, NZ, NXB, NZB, NROT, NZPROF )
      CALL CLZPRF
     I ( DUADV,   1, 12, 1, DTBGT, ZPROF,
     D       NX, NZ, NXB, NZB,    1, NZPROF )
      CALL CLZPRF
     I ( DVADV,   1, 22, 1, DTBGT, ZPROF,
     D       NX, NZ, NXB, NZB,    1, NZPROF )
      CALL CLZPRF
     I ( DUVIS,   1, 13, 1, DTBGT, ZPROF,
     D       NX, NZ, NXB, NZB,    1, NZPROF )
      CALL CLZPRF
     I ( DVVIS,   1, 23, 1, DTBGT, ZPROF,
     D       NX, NZ, NXB, NZB,    1, NZPROF )
      CALL CLZPRF
     I ( DUNLV,   1, 14, 1, DTBGT, ZPROF,
     D       NX, NZ, NXB, NZB,    1, NZPROF )
      CALL CLZPRF
     I ( DVNLV,   1, 24, 1, DTBGT, ZPROF,
     D       NX, NZ, NXB, NZB,    1, NZPROF )
      CALL CLZPRF
     I ( DUCORI,   1, 15, 1, DTBGT, ZPROF,
     D       NX, NZ, NXB, NZB,    1, NZPROF )
      CALL CLZPRF
     I ( DVCORI,   1, 25, 1, DTBGT, ZPROF,
     D       NX, NZ, NXB, NZB,    1, NZPROF )
C     CALL CLZPRF
C    I ( DUFRIC,   1, 16, 1, DTBGT, ZPROF,
C    D       NX, NZ, NXB, NZB,    1, NZPROF )
C     CALL CLZPRF
C    I ( DVFRIC,   1, 26, 1, DTBGT, ZPROF,
C    D       NX, NZ, NXB, NZB,    1, NZPROF )
C
C=====================================================================
C     Kinetic Energy Budget (2001/01/16)

      call EMONIT
     I  (  U    , V    , W    , 
     I     DUVIS, DVVIS, DWVIS,
     I     DUNLV, DVNLV, DWNLV,
     I     DUFRIC,DVFRIC,
     I     BUOY , 
     I     IT1  , IT2  ,
     O     EKVIS, EKNLV, EKBUOY,
     B     DENS0, FDZ  , FDZM  , 
     G     DX   ,  
     D     NX   , NXB  , NZ    , NZB, NROT )

      CALL CLZPRF
     I  (  EKVIS,   1, 31, 1, DTBGT, ZPROF,
     D     NX, NZ, NXB, NZB,    1, NZPROF )
      CALL CLZPRF
     I  (  EKNLV,   1, 32, 1, DTBGT, ZPROF,
     D     NX, NZ, NXB, NZB,    1, NZPROF )
      CALL CLZPRF
     I  (  EKBUOY,  1, 33, 1, DTBGT, ZPROF,
     D     NX, NZ, NXB, NZB,    1, NZPROF )
C
C=====================================================================
C
      CALL CLCONM
     I   ( U, W, IC1, DTI,
     O     CONVM,
     B     DENS0, FDZ,
     D     NX, NZ, NXB, NZB, NROT )
C
C     CALL ADUMP ( CONVM, 'CONVM',  NXB, NZB, 1, 1 , IFDUMP,TIME )
C
      CALL CLPRES
     M   (  DUTOTL(-2,-2), DWTOTL(-2,-2),
     I      CONVM(-2,-2),
     O      PFUNC(-2,-2),  PPAI(-2,-2,IC2),
     W      PDIV (-2,-2),  PDIVM(-2,-2), PFUNCM(-2,-2),
     W      ALPHA(-2,-2), GAMMA(-2,-2),
     M      VMODE(-2,1), TMODE(-2,1), FMODE(-2,1),
     M      TRIGS(1),  IFAC(1), UUFT, ANFT,
     G         DX,
     C      CPDRY,
     B      FDZM, TPOT0,
     D      NX, NZ, NXB, NZB           )
C
      CALL ADUMP ( DUTOTL , 'DU/DT ', NXB, NZB, 1, 1 , IFDUMP,TIME )
C     CALL ADUMP ( DVTOTL , 'DV/DT ', NXB, NZB, 1, 1 , IFDUMP,TIME )
      CALL ADUMP ( DWTOTL , 'DW/DT ', NXB, NZB, 1, 1 , IFDUMP,TIME )
C
C     WRITE(6,*) 'CLPRES END'
C
C---------------------------------------------------------------------
C      Single step time integration.
C---------------------------------------------------------------------
C
      DO 31 IZ = 0, NZ
      DO 31 IX = 1, NX
         U(IX,IZ,IC2) = U(IX,IZ,IC1)  + DTI * DUTOTL(IX,IZ)
         W(IX,IZ,IC2) = W(IX,IZ,IC1)  + DTI * DWTOTL(IX,IZ)
         V(IX,IZ,IC2) = V(IX,IZ,IC1)  + DTI * DVTOTL(IX,IZ)
   31 CONTINUE
C
      CALL BOUND ( U, NX, NZ, NXB, NZB, NROT, 1, 1, IC2 )
      CALL BOUND ( V, NX, NZ, NXB, NZB, NROT, 1, 1, IC2 )
      CALL BOUND ( W, NX, NZ, NXB, NZB, NROT,-1, 0, IC2 )
C
      CALL ADUMP ( U , 'U', NXB, NZB, 2, IC2 , IFDUMP,TIME )
C     CALL ADUMP ( V , 'V', NXB, NZB, 2, IC2 , IFDUMP,TIME )
      CALL ADUMP ( W , 'W', NXB, NZB, 2, IC2 , IFDUMP,TIME )
      CALL ADUMP ( PPAI , 'PPAI', NXB, NZB, 2, IC2 , IFDUMP,TIME )
C
C---------------------------------------------------------------------
C
C     CALL CLCONM
C    I   ( U, W, IC2, DTI,
C    O     CONVM,
C    B     DENS0, FDZ,
C    D     NX, NZ, NXB, NZB, NROT )
C
C     CALL ADUMP ( CONVM,' CONVM NEW', NXB, NZB, 1, 1 )
C
C=====================================================================
C     End of the calculation of the Equations of Motion.
C=====================================================================
C
C
C----------------------------------------------------------------------
C     Turbulent diffusivity and viscosity .
C----------------------------------------------------------------------
C
      CALL CLWFLX
     I   ( U, W,
     I     IT1,
     O     UFLX, WFLX,
     B     FDZ, DENS0,
     D     NX, NZ, NXB, NZB, NROT )
C
C---------------------------------------------------------------------
C
      CALL ADDIF1
     I   (   CDTURB,
     I       UFLX, WFLX,
     I       CD1 ,
     I       IT1, IT2,
     W       FLQX , FLQZ ,
     O       DTQADV, DTQDIF,
     B       DENS0,
     G       DX, FDZ, FDZM,
     D       NX, NZ, NXB, NZB, NROT )
C
      SCTURB = 2000.0D0
      szturb = 0.0d0
C
      CALL NLDIFV
     I   (   CDTURB, SCTURB, szturb, dti, dx,
     I       IT2,
     O       DTQDNL,
     W       QLAP, CDQNLX, CDQNLZ,
     B       zgrid, DENS0, FDZ,
     D       NX, NZ, NXB, NZB, NROT )
C
      CALL CLTURB
     O  (  SSTURB, DTPVIS,
     I     U     , W     , TPOT  , QVAP  , QCLW  , QRAI  , CDTURB,
     I     IT1   , IT2   ,
     B     TPOT0 , QVAP0 , PPAI0 ,
     G     ZGRID , FDZ   , FDZM  , DX    ,
     C     GRAV  , CPDRY , ALATNT, RMVAP ,
     W     STABZ , EDFM1 , EDFM2 ,
     D     NX    , NZ    , NXB   , NZB   , NROT  )
C
C
      DO 41 IZ = 0, NZ-1
      DO 41 IX = 1, NX
         CDTURB(IX,IZ,IC2) =
     \     MAX ( 0.0D0, MIN ( 800.0D0,
     \     ( CDTURB(IX,IZ,IC1)
     \     + DTI * (   DTQADV(IX,IZ) + DTQDIF(IX,IZ) + DTQDNL(IX,IZ)
     \               + SSTURB(IX,IZ) )  ) )  )
   41 CONTINUE
C
C
      CALL BOUND ( CDTURB, NX, NZ, NXB, NZB, NROT, 1, 1, IC2 )
C
      CALL ADUMP ( CDTURB , 'CDTURBL', NXB, NZB, 2, IC2, IFDUMP,TIME )
C=====================================================================
C=====================================================================
C     Thermodynamic equation and
C     Conservation of Water substances.
C---------------------------------------------------------------------
C     Potential temperature.
C---------------------------------------------------------------------
C
      CALL ADDIF0
     I   (   TPOT0(-2),
     I       UFLX(-2,-2), WFLX(-2,-2),
     I       CD1(-2,-2),
     W       QADVM(-2),
     O       DTQAD0(-2,-2), DTQDI0(-2,-2),
     B       DENS0(-2),
     G       DX, FDZ(-2), FDZM(-2),
     D       NX, NZ, NXB, NZB  )
C
      CALL ADDIF1
     I   (   TPOT,
     I       UFLX, WFLX,
     I       CD1 ,
     I       IT1, IT2,
     W       FLQX , FLQZ ,
     O       DTQADV, DTQDIF,
     B       DENS0,
     G       DX, FDZ, FDZM,
     D       NX, NZ, NXB, NZB, NROT )
C
C      SCTPOT = 1.0D0
C      sztpot = 1.0d0 / 10000.0d0
C
c      CALL NLDIFV
c     I   (   TPOT, SCTPOT, sztpot, dti, dx,
c     I       IT2,
c     O       DTQDNL,
c     W       QLAP, CDQNLX, CDQNLZ,
c     B       zgrid, DENS0, FDZ,
c     D       NX, NZ, NXB, NZB, NROT )

      DO IZ = -2, NZB
         DO IX = -2, NXB
            DTQDNL(IX,IZ) = 0.0D0
         END DO
      END DO

C
      CALL ADUMP ( DTQADV , 'DTQADV', NXB, NZB, 1, 1 , IFDUMP,TIME )
      CALL ADUMP ( DTQDIF , 'DTQDIF', NXB, NZB, 1, 1 , IFDUMP,TIME )
      CALL ADUMP ( DTQDNL , 'DTQDNL', NXB, NZB, 1, 1 , IFDUMP,TIME )
      CALL ADUMP ( DTQDI0 , 'DTQDI0', NXB, NZB, 1, 1 , IFDUMP,TIME )
      CALL ADUMP ( DTQAD0 , 'DTQAD0', NXB, NZB, 1, 1 , IFDUMP,TIME )
C
C      write(*,*) 'call CLRAD'

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

      if ( TIME .eq. TIME0 .or. MOD( TIME, 60.0D0 ) .eq. 0.0D0 ) then
C
C *** CO2 Infrared Heaating/Cooling ***
C
      CALL CLRAD
     I   ( TPOT  , TAUN  , PSFC  ,
     I     IT2 ,
     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     IT2  , 
     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
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  , IT2,
     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  , IT2,
     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
            DTQIRAD(ix,iz) = QRAD1   (ix,iz) + QRAD2   (ix,iz)
            FXQRADU(ix,iz) = FXQRADU1(ix,iz) + FXQRADU2(ix,iz)
            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 iz = -2, NZB
         do ix = -2, NXB
            DTPRAD(ix,iz) = DTIRAD(ix,iz) + DTSRAD(ix,iz) +
     \           DTQIRAD(ix,iz) + DTQSRAD(ix,iz)
         end do
      end do

      end if

C      DO 51 IZ = 0, NZ-1
      DO 51 IZ = 0, NZ
      DO 51 IX = 1, NX
         TPOT1(IX,IZ) = TPOT(IX,IZ,IC1)
     \     + DTI  * (   DTQADV(IX,IZ) + DTQDIF(IX,IZ) + DTQDNL(IX,IZ)
     \                + DTQAD0(IX,IZ) + DTQDI0(IX,IZ)
     \                + DTPRAD(IX,IZ) + DTPVIS(IX,IZ)  )
   51 CONTINUE
      
! N 層目は放射計算だけ別途行なう

C      DO 510 IX = 1, NX
C         TPOT(IX,NZ,IC2) = TPOT(IX,NZ,IC1) + DTI * DTPRAD(IX,NZ) 
C  510 CONTINUE

      DO 511 IX = 1, NX
         TPOT1(IX,0) = TPOT1(IX,0)
     \     + DTI * FHEAT(IX,1) / DX / FDZ(0) / DENS0(0) / CPDRY
         TPOT1(IX,NZ-1) = TPOT1(IX,NZ-1)
     \     + DTI * FHEAT(IX,2) / DX / FDZ(NZ-1) / DENS0(NZ-1) / CPDRY
  511 CONTINUE
C
C
      DO 501 IX = 1, NX
c         TPOT1(IX,NZ-1) = 0.0D0
  501 CONTINUE
C
      CALL BOUND ( TPOT1, NX, NZ, NXB, NZB, 1, 1, 1, 1 )
C
      CALL CLZPRF
     I ( TPOT , IC1, 41, 0, 0.0D0, ZPROF,
     D       NX, NZ, NXB, NZB, NROT, NZPROF )
      CALL CLZPRF
     I ( DTQADV,   1, 42, 1, DTBGT, ZPROF,
     D       NX, NZ, NXB, NZB,    1, NZPROF )
      CALL CLZPRF
     I ( DTQDIF,   1, 43, 1, DTBGT, ZPROF,
     D       NX, NZ, NXB, NZB,    1, NZPROF )
      CALL CLZPRF
     I ( DTQDNL,   1, 44, 1, DTBGT, ZPROF,
     D       NX, NZ, NXB, NZB,    1, NZPROF )
      CALL CLZPRF
     I ( DTQAD0,   1, 45, 1, DTBGT, ZPROF,
     D       NX, NZ, NXB, NZB,    1, NZPROF )
      CALL CLZPRF
     I ( DTQDI0,   1, 46, 1, DTBGT, ZPROF,
     D       NX, NZ, NXB, NZB,    1, NZPROF )
      CALL CLZPRF
     I ( DTPVIS,   1, 47, 1, DTBGT, ZPROF,
     D       NX, NZ, NXB, NZB,    1, NZPROF )
      CALL CLZPRF
     I ( DTIRAD,   1, 48, 1, DTBGT, ZPROF,
     D       NX, NZ, NXB, NZB,    1, NZPROF )
      CALL CLZPRF
     I ( DTSRAD,   1, 49, 1, DTBGT, ZPROF,
     D       NX, NZ, NXB, NZB,    1, NZPROF )
C     CALL CLZPRF
C    I ( FHEAT ,   1, 48, 1, DTBGT, ZPROF,
C    D       NX, NZ, NXB, NZB,    1, NZPROF )
C
C
C---------------------------------------------------------------------
C     Surface Temp.
C---------------------------------------------------------------------

      do ix = -2, NXB
C         F0(ix) = S0 * COSZ - ( SXRADD(ix,NZ+1) - SXRADD(ix,0) )
         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      if ( TIME .eq. TIME0 .or. MOD( TIME, 60.0D0 ) .eq. 0.0D0 ) then

      call CLTSFC
     I   ( FHEAT , FXRADG, F0,
     O     TEMPG , TSFC  ,
     C     SIGMA , ALBEDO, CPGRND, 
     C     CDTURBG,DTI   , DENSG ,
     B     GDZ   , GDZM  , 
     D     NX    , NXB   , NZ   , NZB   ,
     D     NZG   , NZGB                   )

C      end if
         
C---------------------------------------------------------------------
C     Water Vapor.
C---------------------------------------------------------------------
C
      CALL ADDIF0
     I   (   QVAP0(-2),
     I       UFLX(-2,-2), WFLX(-2,-2),
     I       CD1(-2,-2),
     W       QADVM(-2),
     O       DTQAD0(-2,-2), DTQDI0(-2,-2),
     B       DENS0(-2),
     G       DX, FDZ(-2), FDZM(-2),
     D       NX, NZ, NXB, NZB  )
C
      CALL ADDIF1
     I   (   QVAP,
     I       UFLX, WFLX,
     I       CD1 ,
     I       IT1, IT2,
     W       FLQX , FLQZ ,
     O       DTQADV, DTQDIF,
     B       DENS0,
     G       DX, FDZ, FDZM,
     D       NX, NZ, NXB, NZB, NROT )

      call ADFAL1
     I   (   QVAP, IT1, 
     W       WFLX,
     O       DTQFAL, 
     B       DENS0, PRES0, 
     C       GRAV , QDENS, QRMOD, 
     G       DX, FDZ, FDZM,
     D       NX, NZ, NXB, NZB, NROT )
C
      SCQVAP =  2.0D-3
      szqvap = - 1.0d0 / 3000.0d0
C
      CALL NLDIFV
     I   (   QVAP, SCQVAP, szqvap, dti, dx,
     I       IT2,
     O       DTQDNL,
     W       QLAP, CDQNLX, CDQNLZ,
     B       zgrid, DENS0, FDZ,
     D       NX, NZ, NXB, NZB, NROT )
C
      CALL ADUMP ( DTQADV , 'DQVADV', NXB, NZB, 1, 1 , IFDUMP,TIME )
      CALL ADUMP ( DTQDIF , 'DQVDIF', NXB, NZB, 1, 1 , IFDUMP,TIME )
      CALL ADUMP ( DTQDNL , 'DQVDNL', NXB, NZB, 1, 1 , IFDUMP,TIME )
C     CALL ADUMP ( DTQDI0 , 'DQVDI0', NXB, NZB, 1, 1 , IFDUMP,TIME )
C     CALL ADUMP ( DTQAD0 , 'DQVAD0', NXB, NZB, 1, 1 , IFDUMP,TIME )
C
C
      DO 52 IZ = 0, NZ-1
      DO 52 IX = 1, NX
         QVAP1(IX,IZ) = QVAP(IX,IZ,IC1)
     \     + DTI  * (   DTQADV(IX,IZ) + DTQDIF(IX,IZ) + DTQDNL(IX,IZ)
     \                + DTQAD0(IX,IZ) + DTQDI0(IX,IZ) + DTQFAL(IX,IZ) )
   52 CONTINUE
C
      DO 521 IX = 1, NX
         QVAP1(IX,0) = QVAP1(IX,0)
     \     + DTI * FEVAP(IX,1) / DX / FDZ(0) / DENS0(0)
         QVAP1(IX,NZ-1) = QVAP1(IX,NZ-1)
     \     + DTI * FEVAP(IX,2) / DX / FDZ(NZ-1) / DENS0(NZ-1)
  521 CONTINUE
C
      CALL BOUND ( QVAP1, NX, NZ, NXB, NZB, 1, 1, 1, 1 )
C
      CALL QNFILL
     I   ( QVAP1, 1,
     O     DQFILL,
     W     QSUMPN,
     B     DENS0, FDZ,
     D     NX, NZ, NXB, NZB, 1 )
C
      CALL ADUMP ( DQFILL , 'DQVFIL', NXB, NZB, 1, 1 , IFDUMP,TIME )
C
      CALL CLZPRF
     I ( QVAP , IC1, 51, 0, 0.0D0, ZPROF,
     D       NX, NZ, NXB, NZB, NROT, NZPROF )
      CALL CLZPRF
     I ( DTQADV,   1, 52, 1, DTBGT, ZPROF,
     D       NX, NZ, NXB, NZB,    1, NZPROF )
      CALL CLZPRF
     I ( DTQDIF,   1, 53, 1, DTBGT, ZPROF,
     D       NX, NZ, NXB, NZB,    1, NZPROF )
      CALL CLZPRF
     I ( DTQDNL,   1, 54, 1, DTBGT, ZPROF,
     D       NX, NZ, NXB, NZB,    1, NZPROF )
      CALL CLZPRF
     I ( DTQAD0,   1, 55, 1, DTBGT, ZPROF,
     D       NX, NZ, NXB, NZB,    1, NZPROF )
      CALL CLZPRF
     I ( DTQDI0,   1, 56, 1, DTBGT, ZPROF,
     D       NX, NZ, NXB, NZB,    1, NZPROF )
      CALL CLZPRF
     I ( DQFILL,   1, 57, 1, DTBGT/DTI, ZPROF,
     D       NX, NZ, NXB, NZB,    1, NZPROF )
      CALL CLZPRF
     I ( DTQFAL,   1, 58, 1, DTBGT/DTI, ZPROF,
     D       NX, NZ, NXB, NZB,    1, NZPROF )
C     CALL CLZPRF
C    I ( FEVAP ,   1, 59, 1, DTBGT, ZPROF,
C    D       NX, NZ, NXB, NZB,    1, NZPROF )
C
C----------------------------------------------------------------------
C     Cloud Water.
C----------------------------------------------------------------------
C
      CALL ADUMP ( QCLW, 'QCLW 1', NXB, NZB, NROT, IT1 , IFDUMP,TIME )
      CALL ADUMP ( QCLW, 'QCLW 2', NXB, NZB, NROT, IT2 , IFDUMP,TIME )

      CALL ADDIF1
     I   (   QCLW,
     I       UFLX, WFLX,
     I       CD1 ,
     I       IT1, IT2,
     W       FLQX , FLQZ ,
     O       DTQADV, DTQDIF,
     B       DENS0,
     G       DX, FDZ, FDZM,
     D       NX, NZ, NXB, NZB, NROT )
C
      SCQCLW = 0.5D-3
      szqclw = -1.0d0 / 3000.0d0
C
      CALL NLDIFV
     I   (   QCLW, SCQCLW, szqclw, dti, dx, 
     I       IT2,
     O       DTQDNL,
     W       QLAP, CDQNLX, CDQNLZ,
     B       zgrid, DENS0, FDZ,
     D       NX, NZ, NXB, NZB, NROT )
C
      CALL ADUMP ( DTQADV , 'DQCADV', NXB, NZB, 1, 1 , IFDUMP,TIME )
      CALL ADUMP ( DTQDIF , 'DQCDIF', NXB, NZB, 1, 1 , IFDUMP,TIME )
      CALL ADUMP ( DTQDNL , 'DQCDNL', NXB, NZB, 1, 1 , IFDUMP,TIME )
C
C
      DO 53 IZ = 0, NZ-1
      DO 53 IX = 1, NX
         QCLW1(IX,IZ) = QCLW(IX,IZ,IC1)
     \     + DTI  * ( DTQADV(IX,IZ) + DTQDIF(IX,IZ) + DTQDNL(IX,IZ) )
   53 CONTINUE
C
      CALL BOUND ( QCLW1, NX, NZ, NXB, NZB, 1, 1, 1, 1 )
C
C
      CALL QNFILL
     I   ( QCLW1, 1,
     O     DQFILL,
     W     QSUMPN,
     B     DENS0, FDZ,
     D     NX, NZ, NXB, NZB, 1 )
C
      CALL ADUMP ( DQFILL , 'DQCFIL', NXB, NZB, 1, 1 , IFDUMP,TIME )
C
      CALL CLZPRF
     I ( QCLW , IC1, 61, 0, 0.0D0, ZPROF,
     D       NX, NZ, NXB, NZB, NROT, NZPROF )
      CALL CLZPRF
     I ( DTQADV,   1, 62, 1, DTBGT, ZPROF,
     D       NX, NZ, NXB, NZB,    1, NZPROF )
      CALL CLZPRF
     I ( DTQDIF,   1, 63, 1, DTBGT, ZPROF,
     D       NX, NZ, NXB, NZB,    1, NZPROF )
      CALL CLZPRF
     I ( DTQDNL,   1, 64, 1, DTBGT, ZPROF,
     D       NX, NZ, NXB, NZB,    1, NZPROF )
      CALL CLZPRF
     I ( DQFILL,   1, 65, 1, DTBGT/DTI, ZPROF,
     D       NX, NZ, NXB, NZB,    1, NZPROF )
C
C----------------------------------------------------------------------
C     Rain Water.
C----------------------------------------------------------------------
C
      CALL ADDIF1
     I   (   QRAI,
     I       UFLX, WFLX,
     I       CD1 ,
     I       IT1, IT2,
     W       FLQX , FLQZ ,
     O       DTQADV, DTQDIF,
     B       DENS0,
     G       DX, FDZ, FDZM,
     D       NX, NZ, NXB, NZB, NROT )
C
      SCQRAI = 1.0D-3
      szqrai = -1.0d0  / 3000.0d0
C
      CALL NLDIFV
     I   (   QRAI, SCQRAI, szqrai, dti, dx,
     I       IT2,
     O       DTQDNL,
     W       QLAP, CDQNLX, CDQNLZ,
     B       zgrid, DENS0, FDZ,
     D       NX, NZ, NXB, NZB, NROT )
C
      CALL ADUMP ( DTQADV , 'DQRADV', NXB, NZB, 1, 1 , IFDUMP,TIME )
      CALL ADUMP ( DTQDIF , 'DQRDIF', NXB, NZB, 1, 1 , IFDUMP,TIME )
      CALL ADUMP ( DTQDNL , 'DQRDNL', NXB, NZB, 1, 1 , IFDUMP,TIME )
C
C
      DO 54 IZ = 0, NZ-1
      DO 54 IX = 1, NX
         QRAI1(IX,IZ) = QRAI(IX,IZ,IC1)
     \     + DTI  * ( DTQADV(IX,IZ) + DTQDIF(IX,IZ) + DTQDNL(IX,IZ) )
   54 CONTINUE
C
      CALL CLZPRF
     I ( QRAI , IC1, 71, 0, 0.0D0, ZPROF,
     D       NX, NZ, NXB, NZB, NROT, NZPROF )
      CALL CLZPRF
     I ( DTQADV,   1, 72, 1, DTBGT, ZPROF,
     D       NX, NZ, NXB, NZB,    1, NZPROF )
      CALL CLZPRF
     I ( DTQDIF,   1, 73, 1, DTBGT, ZPROF,
     D       NX, NZ, NXB, NZB,    1, NZPROF )
      CALL CLZPRF
     I ( DTQDNL,   1, 74, 1, DTBGT, ZPROF,
     D       NX, NZ, NXB, NZB,    1, NZPROF )
C
      CALL FARAIN2
     I   (  QRAI ,
     I      DTI  , IT2 ,
     M      QRAI1,
     W      DQFILL, QSUMPN,
     O      VRFALL, FXRAIN, FWPREC,
     B      DENS0,
     G      DX   , FDZ  ,
     D      NX, NZ, NXB, NZB, NROT )
C
      CALL ADUMP ( FXRAIN , 'FXRAIN', NXB, NZB, 1, 1 , IFDUMP,TIME )
      CALL ADUMP ( DQFILL , 'DQRFIL', NXB, NZB, 1, 1 , IFDUMP,TIME )
C
      CALL CLZPRF
     I ( FXRAIN , 1, 91, 1, DTBGT, ZPROF,
     D       NX, NZ, NXB, NZB, 1, NZPROF )
C
      CALL CLZPRF
     I ( DQFILL,   1, 75, 1, DTBGT/DTI, ZPROF,
     D       NX, NZ, NXB, NZB,    1, NZPROF )
C
C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C      Loss of liquid water due to precipitation to ground.
C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
      DO 6000 IX = 1, NX
         WSLIQD = WSLIQD - FWPREC(IX,1) * DTBGT/DTI * DX
     &                   - FWPREC(IX,2) * DTBGT/DTI * DX
         FWPREA(IX,1) = FWPREA(IX,1) + FWPREC(IX,1) * DTBGT/DTI
         FWPREA(IX,2) = FWPREA(IX,2) + FWPREC(IX,2) * DTBGT/DTI
 6000 CONTINUE
C
C +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C      Dissipation of the Potential snergy of water substance
C          via precipitation of Rainwater.
C +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
      DO 7000 IZ = 1, NZ-1
         DMZ = DTI*FDZ(IZ)*(DENS0(IZ)+DENS0(IZ-1))/2.0D0 * DX**2
      DO 7000 IX = 1, NX
         DPWFAL = DPWFAL + DMZ * FXRAIN(IX,IZ)
 7000 CONTINUE
C
C +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
      CALL CLPHYS
     I   (  TPOT1, QVAP1, QCLW1, QRAI1,
     I      DTI  , IT2  , IC2  ,
     M      TPOT , QVAP , QCLW , QRAI ,
     O      PCCON, PVEVC, PVEVR,
     O      PRMIC, PCONV, PCOAG,
     B      DENS0, PPAI0, QVAP0, pres0, tpot0,
     C      CPDRY, ALATNT, rmvap,
     D      NX, NZ, NXB, NZB, NROT  )
C
      CALL ADUMP ( TPOT1  , 'TPOT1 ', NXB, NZB, 1, 1 , IFDUMP,TIME )
      CALL ADUMP ( QVAP1  , 'QVAP1 ', NXB, NZB, 1, 1 , IFDUMP,TIME )
      CALL ADUMP ( QCLW1  , 'QCLW1 ', NXB, NZB, 1, 1 , IFDUMP,TIME )
      CALL ADUMP ( QRAI1  , 'QRAI1 ', NXB, NZB, 1, 1 , IFDUMP,TIME )
      CALL ADUMP ( PCCON  , 'PCCON ', NXB, NZB, 1, 1 , IFDUMP,TIME )
      CALL ADUMP ( PVEVC  , 'PVEVC ', NXB, NZB, 1, 1 , IFDUMP,TIME )
      CALL ADUMP ( PVEVR  , 'PVEVR ', NXB, NZB, 1, 1 , IFDUMP,TIME )
      CALL ADUMP ( PRMIC  , 'PRMIC ', NXB, NZB, 1, 1 , IFDUMP,TIME )
      CALL ADUMP ( PCONV  , 'PCONV ', NXB, NZB, 1, 1 , IFDUMP,TIME )
      CALL ADUMP ( PCOAG  , 'PCOAG ', NXB, NZB, 1, 1 , IFDUMP,TIME )
C
      CALL CLZPRF
     I ( PCCON ,   1, 92, 1, DTBGT/DTI, ZPROF,
     D       NX, NZ, NXB, NZB,    1, NZPROF )
      CALL CLZPRF
     I ( PVEVC ,   1, 93, 1, DTBGT/DTI, ZPROF,
     D       NX, NZ, NXB, NZB,    1, NZPROF )
      CALL CLZPRF
     I ( PVEVR ,   1, 94, 1, DTBGT/DTI, ZPROF,
     D       NX, NZ, NXB, NZB,    1, NZPROF )
      CALL CLZPRF
     I ( PCONV ,   1, 95, 1, DTBGT/DTI, ZPROF,
     D       NX, NZ, NXB, NZB,    1, NZPROF )
      CALL CLZPRF
     I ( PCOAG ,   1, 96, 1, DTBGT/DTI, ZPROF,
     D       NX, NZ, NXB, NZB,    1, NZPROF )
C +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C     Conversion among water substances.
C +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
      DO 8000 IZ = 0, NZ-1
         DMZ = FDZ(IZ) * DENS0(IZ) * DX**2 * DTBGT / DTI
      DO 8000 IX = 1, NX
         WCVTOL = WCVTOL +
     \      DMZ * ( PCCON(IX,IZ) - PVEVC(IX,IZ) - PVEVR(IX,IZ) )
 8000 CONTINUE
C
C +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C      Conversion from Latent to Thermal energy
C +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
      DO 9000 IZ = 0, NZ-1
         DMZ = FDZ(IZ) * DENS0(IZ) * DX**2  * ALATNT * DTBGT / DTI
      DO 9000 IX = 1, NX
         CLATHR = CLATHR +
     \      DMZ * ( PCCON(IX,IZ) - PVEVC(IX,IZ) - PVEVR(IX,IZ) )
 9000 CONTINUE
C
C +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C      Generation of the Potential energy of Water substance.
C +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
      ZREF = 3.0D+4
C
      DO 9010 IZ = 0, NZ-1
         DMZL =   FDZ(IZ) * DENS0(IZ) * DX**2  * ZGRID(IZ) * GRAV
     \          * DTBGT / DTI
         DMZV = - FDZ(IZ) * DENS0(IZ) * DX**2  * ( ZGRID(IZ) - ZREF )
     \                    * GRAV * ( 1.0D0/RMVAP - 1.0D0 )
     \                    * DTBGT / DTI
      DO 9010 IX = 1, NX
         SEPOTW = SEPOTW +
     \      DMZL * (   PCCON(IX,IZ) - PVEVC(IX,IZ) - PVEVR(IX,IZ) )
     \    + DMZV * ( - PCCON(IX,IZ) + PVEVC(IX,IZ) + PVEVR(IX,IZ) )
 9010 CONTINUE
C
      DO 9020 IX = 1, NX
         SEPOTW = SEPOTW
     \      + ZREF * GRAV * ( 1.0D0/RMVAP - 1.0D0 )
     \             * FEVAP(IX,1) * DTBGT
     \      + (ZREF-ZGRID(NZ-1)) * GRAV * ( 1.0D0/RMVAP - 1.0D0 )
     \             * FEVAP(IX,2) * DTBGT
 9020 CONTINUE
C
C---------------------------------------------------------------------
C
      DO 8800 IZ = 0, NZ-1
      DO 8800 IX = 1, NX
         XQFILL = XQFILL
     \     - FDZ(IZ) * DENS0(IZ) * DX**2 * DTBGT / DTI
     \               * (   MIN( 0.0D0, QVAP(IX,IZ,IC2) )
     \                   + MIN( 0.0D0, QCLW(IX,IZ,IC2) )
     \                   + MIN( 0.0D0, QRAI(IX,IZ,IC2) ) )
 8800 CONTINUE
C
C
      DO 8801 IZ = 0, NZ-1
      DO 8801 IX = 1, NX
         QRAI(IX,IZ,IC2) = MAX( 0.0D0, QRAI(IX,IZ,IC2) )
         QCLW(IX,IZ,IC2) = MAX( 0.0D0, QCLW(IX,IZ,IC2) )
         QVAP(IX,IZ,IC2) = MAX( 0.0D0, QVAP(IX,IZ,IC2) )
 8801 CONTINUE
C
      CALL BOUND ( QRAI, NX, NZ, NXB, NZB, NROT, 1, 1, IC2 )
      CALL BOUND ( QCLW, NX, NZ, NXB, NZB, NROT, 1, 1, IC2 )
      CALL BOUND ( QVAP, NX, NZ, NXB, NZB, NROT, 1, 1, IC2 )
C
C=====================================================================
C
 9876 RETURN
      END
