C**********************************************************************
C
      PROGRAM MOIST
C
C**********************************************************************
C
      IMPLICIT REAL*8 ( A-H, O-Z )
C
C**********************************************************************
C
      include 'grid_size.f'
C
      PARAMETER ( NZPROF = 100 )
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 )
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 *** 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
      NAMELIST /CONST2/ CDRAG, VSFC0, PSFC, TSFC, QVSFC, FCOLI
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 = 41, 48 )
     \   / 'TPMEAN', 'DTPADV', 'DTPDIF', 'DTPDNL', 'DTPAD0', 'DTPDI0',
     \     'DTPRAD', 'FHEAT' /
      DATA  ( NMZPRF(IST), IST = 51, 58 )
     \   / 'QVMEAN', 'DQVADV', 'DQVDIF', 'DQVDNL', 'DQVAD0', 'DQVDI0',
     \     'DQVNFL', '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
      WRITE(16, MODEL  )
      WRITE(16, CONST1 )
      WRITE(16, CONST2 )
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
      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
      WRITE(6,*)'TIME0, NLOOP1, NLOOP2, DTIME, NTMONI, NTKUBU'
C
C
      READ (5,*) TIME0, NLOOP1, NLOOP2, DTIME, NTMONI , NTKUBU
C
C
      WRITE(6,*) TIME0, NLOOP1, NLOOP2, DTIME, NTMONI , NTKUBU
      WRITE(16,*) TIME0, NLOOP1, NLOOP2, DTIME, NTMONI , NTKUBU
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
      ELSE
C
      CALL FLREAD
     I   (  TIME0,  1,
     O      TIME,
     O      U, V, W, TPOT, QVAP, QCLW, QRAI, PPAI, CDTURB,
     D      NX, NZ, NXB, NZB, NROT           )
C
      END IF

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,
     D         NX, NZ, NXB, NZB, 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,
     C     DX, GRAV, CPDRY, ALATNT, RMVAP,
     C     FCOLI, CDRAG, VSFC0, TSFC, QVSFC,
     V     U   , V   , W   , PPAI, TPOT, QVAP, QCLW, QRAI, CDTURB,
     F     FUMOM, FVMOM, FHEAT, FEVAP, FWPREC,
     F     FUMOMA, FVMOMA, FHEATA, FEVAPA, FWPREA,
     B     ZGRID, FDZ  , FDZM , TEMP0, TPOT0,
     B     DENS0, PPAI0, QVAP0, pres0,
     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,
     C     DX, GRAV, CPDRY, ALATNT, RMVAP,
     C     FCOLI, CDRAG, VSFC0, TSFC, QVSFC,
     V     U   , V   , W   , PPAI, TPOT, QVAP, QCLW, QRAI, CDTURB,
     F     FUMOM, FVMOM, FHEAT, FEVAP, FWPREC,
     F     FUMOMA, FVMOMA, FHEATA, FEVAPA, FWPREA,
     B     ZGRID, FDZ  , FDZM , TEMP0, TPOT0,
     B     DENS0, PPAI0, QVAP0, pres0,
     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
      WRITE(61) TIME, NMZPRF, 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,
     D      NX, NZ, NXB, NZB, NROT           )
C
C
 9999 CONTINUE
C
C
      CALL FCLOSE
C
 9091 CONTINUE
C
      STOP
      END
      SUBROUTINE TSTEP
     I   ( IT1, IT2, IC1, IC2, DTI,
     I     TIME, IFDUMP,
     C     DX, GRAV, CPDRY, ALATNT, RMVAP,
     C     FCOLI, CDRAG, VSFC0, TSFC, QVSFC,
     V     U   , V   , W   , PPAI, TPOT, QVAP, QCLW, QRAI, CDTURB,
     F     FUMOM, FVMOM, FHEAT, FEVAP, FWPREC,
     F     FUMOMA, FVMOMA, FHEATA, FEVAPA, FWPREA,
     B     ZGRID, FDZ  , FDZM , TEMP0, TPOT0,
     B     DENS0, PPAI0, QVAP0, pres0,
     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.f'
C
      PARAMETER ( NZPROF = 100 )
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 )
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 )
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 )
C
      DIMENSION
     \     DQFILL ( -2:NXB, -2:NZB ), QSUMPN ( -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 *** 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( -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, 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
      CALL FXSURF2
     I   ( U,     V,     TPOT , QVAP ,
     I     IT2  ,
     O     FUMOM, FVMOM, FHEAT, FEVAP,
     B     DENS0, TPOT0, QVAP0, PPAI0,
     C     CPDRY,
     S     CDRAG, 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
      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----------------------------------------------------------------------
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,
     I       IT2,
     O       DTQDNL,
     W       QLAP, CDQNLX, CDQNLZ,
     B       zgrid, DENS0, FDZ,
     D       NX, NZ, NXB, NZB, NROT )
C
      CALL CLTURB
     O  (  SSTURB,
     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
      SCTPOT = 1.0D0
      sztpot = 1.0d0 / 10000.0d0
C
      CALL NLDIFV
     I   (   TPOT, SCTPOT, sztpot,
     I       IT2,
     O       DTQDNL,
     W       QLAP, CDQNLX, CDQNLZ,
     B       zgrid, DENS0, FDZ,
     D       NX, NZ, NXB, NZB, NROT )
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
      CALL CLRAD
     I   ( TPOT,
     I     IT2,
     O     DTPRAD,
     W     TPMEAN, ALMEAN, ALDEVI, COOLZ,
     B     TPOT0, TEMP0, ZGRID,
     D     NX, NZ, NXB, NZB, NROT  )
C
      DO 51 IZ = 0, NZ-1
      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)                  )
   51 CONTINUE
C
      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 ( DTPRAD,   1, 47, 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     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 )
C
      SCQVAP =  2.0D-3
      szqvap = - 1.0d0 / 3000.0d0
C
      CALL NLDIFV
     I   (   QVAP, SCQVAP, szqvap,
     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)                )
   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 )
C     CALL CLZPRF
C    I ( FEVAP ,   1, 58, 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,
     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,
     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
