!        Generated by TAPENADE     (INRIA, Tropics team)
!  Tapenade 3.7 (r4786) - 21 Feb 2013 15:53
!
!  Differentiation of ducu in forward (tangent) mode (with options r8):
!   variations   of useful results: raincv rthcuten pratec rqvcuten
!   with respect to varying inputs: th raincv t rthcuten pcps qv
!                z pratec rqvcuten rho dz8w
MODULE g_module_cu_du
   USE module_wrf_error
  REAL    , PARAMETER :: cincap = -10.
  REAL    , PARAMETER :: capemin = 10.
  REAL    , PARAMETER :: dpthmin = 1000.
  REAL    , PARAMETER :: alpha = 0.00002
  REAL    , PARAMETER :: eps = 0.5
  REAL    , PARAMETER :: Vfall = 5.
!--------------------------------------------------------------------
CONTAINS

SUBROUTINE DUCU_D(ids, ide, jds, jde, kds, kde, ims, ime, jms, jme, kms&
&  , kme, its, ite, jts, jte, kts, kte, dt, ktau, dx, rho, rhod, raincv, &
&  raincvd, nca, pratec, pratecd, u, v, th, thd, t, td, w, dz8w, dz8wd, z&
&  , zd, pcps, pcpsd, pi, w0avg, cp, rd, rv, g, xlv, ep2, svp1, svp2, &
&  svp3, svpt0, stepcu, cu_act_flag, warm_rain, cutop, cubot, qv, qvd, &
&  rthcuten, rthcutend, rqvcuten, rqvcutend)
  IMPLICIT NONE
!
!-------------------------------------------------------------
  INTEGER, INTENT(IN) :: ids, ide, jds, jde, kds, kde, ims, ime, jms, &
&  jme, kms, kme, its, ite, jts, jte, kts, kte
  INTEGER, INTENT(IN) :: stepcu
  LOGICAL, INTENT(IN) :: warm_rain
  REAL, INTENT(IN) :: xlv
  REAL, INTENT(IN) :: cp, rd, rv, g, ep2
  REAL, INTENT(IN) :: svp1, svp2, svp3, svpt0
  INTEGER, INTENT(IN) :: ktau
  REAL, DIMENSION(ims:ime, kms:kme, jms:jme), INTENT(IN) :: u, v, w, th&
&  , t, qv, dz8w, z, pcps, rho, pi
  REAL, DIMENSION(ims:ime, kms:kme, jms:jme), INTENT(IN) :: thd, td, qvd&
&  , dz8wd, zd, pcpsd, rhod
!
  REAL, DIMENSION(ims:ime, kms:kme, jms:jme), INTENT(INOUT) :: w0avg
  REAL, INTENT(IN) :: dt, dx
!
  REAL, DIMENSION(ims:ime, jms:jme), INTENT(INOUT) :: raincv, pratec
  REAL, DIMENSION(ims:ime, jms:jme), INTENT(INOUT) :: raincvd, pratecd
  REAL, DIMENSION(ims:ime, jms:jme), INTENT(INOUT) :: nca
  REAL, DIMENSION(ims:ime, jms:jme), INTENT(OUT) :: cubot, cutop
  LOGICAL, DIMENSION(ims:ime, jms:jme), INTENT(INOUT) :: cu_act_flag
!
! Optional arguments
!
  REAL, DIMENSION(ims:ime, kms:kme, jms:jme), OPTIONAL, INTENT(INOUT) ::&
&  rthcuten, rqvcuten
  REAL, DIMENSION(ims:ime, kms:kme, jms:jme), OPTIONAL, INTENT(INOUT) ::&
&  rthcutend, rqvcutend
!
! LOCAL VARS
  LOGICAL :: flag_qr, flag_qi, flag_qs
  REAL, DIMENSION(kts:kte) :: u1d, v1d, t1d, th1d, dz1d, z1d, qv1d, p1d&
&  , rho1d, w0avg1d
  REAL, DIMENSION(kts:kte) :: t1dd, th1dd, dz1dd, z1dd, qv1dd, p1dd, &
&  rho1dd
  REAL, DIMENSION(kts:kte) :: dqvdt, dthdt
  REAL, DIMENSION(kts:kte) :: dqvdtd, dthdtd
  REAL :: pprate, tst, tv, prs, rhoe, w0, scr1, dxsq, rthcumax
  REAL :: pprated
  INTEGER :: i, j, k, i_start, i_end, j_start, j_end, sz, ntst, icldck
!
  dxsq = dx*dx
  ntst = stepcu
  icldck = MOD(ktau, ntst)
  IF (icldck .EQ. 0 .OR. ktau .EQ. 1) THEN
!
!  Keep away from specified and relaxation zone (should be for just specified and nested bc)
    sz = 1
    IF (ids + sz .LT. its) THEN
      i_start = its
    ELSE
      i_start = ids + sz
    END IF
    IF (ide - 1 - sz .GT. ite) THEN
      i_end = ite
    ELSE
      i_end = ide - 1 - sz
    END IF
    IF (jds + sz .LT. jts) THEN
      j_start = jts
    ELSE
      j_start = jds + sz
    END IF
    IF (jde - 1 - sz .GT. jte) THEN
      j_end = jte
    ELSE
      j_end = jde - 1 - sz
    END IF
    qv1dd = 0.0_8
    dthdtd = 0.0_8
    th1dd = 0.0_8
    rho1dd = 0.0_8
    p1dd = 0.0_8
    z1dd = 0.0_8
    dz1dd = 0.0_8
    t1dd = 0.0_8
    dqvdtd = 0.0_8
!
    DO j=j_start,j_end
      DO i=i_start,i_end
        DO k=kts,kte
          dqvdtd(k) = 0.0_8
          dqvdt(k) = 0.
          dthdtd(k) = 0.0_8
          dthdt(k) = 0.
        END DO
        raincvd(i, j) = 0.0_8
        raincv(i, j) = 0.
        pratecd(i, j) = 0.0_8
        pratec(i, j) = 0.
        cutop(i, j) = kts
        cubot(i, j) = kte + 1
!
! assign vars from 3D to 1D
        DO k=kts,kte
          u1d(k) = u(i, k, j)
          v1d(k) = v(i, k, j)
          t1dd(k) = td(i, k, j)
          t1d(k) = t(i, k, j)
          th1dd(k) = thd(i, k, j)
          th1d(k) = th(i, k, j)
          rho1dd(k) = rhod(i, k, j)
          rho1d(k) = rho(i, k, j)
          qv1dd(k) = qvd(i, k, j)
          qv1d(k) = qv(i, k, j)
          IF (qv1d(k) .LT. 1.e-08) THEN
            qv1dd(k) = 0.0_8
            qv1d(k) = 1.e-08
          END IF
          p1dd(k) = pcpsd(i, k, j)
          p1d(k) = pcps(i, k, j)
          w0avg1d(k) = w0avg(i, k, j)
          dz1dd(k) = dz8wd(i, k, j)
          dz1d(k) = dz8w(i, k, j)
          z1dd(k) = zd(i, k, j)
          z1d(k) = z(i, k, j)
        END DO
        CALL DUCU1D_D(i, j, u1d, v1d, t1d, t1dd, qv1d, qv1dd, p1d, p1dd&
&                , dz1d, dz1dd, z1d, z1dd, w0avg1d, dt, dx, dxsq, rho1d, &
&                rho1dd, th1d, th1dd, xlv, cp, rd, rv, g, ep2, svp1, svp2&
&                , svp3, svpt0, dqvdt, dqvdtd, dthdt, dthdtd, pprate, &
&                pprated, nca, ntst, cutop, cubot, ids, ide, jds, jde, &
&                kds, kde, ims, ime, jms, jme, kms, kme, its, ite, jts, &
&                jte, kts, kte)
        IF (PRESENT(rthcuten) .AND. PRESENT(rqvcuten)) THEN
          DO k=kts,kte
            rthcutend(i, k, j) = dthdtd(k)
            rthcuten(i, k, j) = dthdt(k)
            rqvcutend(i, k, j) = dqvdtd(k)
            rqvcuten(i, k, j) = dqvdt(k)
          END DO
          pratecd(i, j) = pprated
          pratec(i, j) = pprate
          raincvd(i, j) = dt*pprated
          raincv(i, j) = pprate*dt
        END IF
      END DO
    END DO
  END IF
END SUBROUTINE DUCU_D

!        Generated by TAPENADE     (INRIA, Tropics team)
!  Tapenade 3.7 (r4786) - 21 Feb 2013 15:53
!
!  Differentiation of ducu1d in forward (tangent) mode (with options r8):
!   variations   of useful results: dthdt pprate dqvdt
!   with respect to varying inputs: dthdt p0 z t0 th0 rhoe dzq
!                dqvdt qv0
! ****************************************************************************
!-----------------------------------------------------------
SUBROUTINE DUCU1D_D(i, j, u0, v0, t0, t0d, qv0, qv0d, p0, p0d, dzq, dzqd&
&  , z, zd, w0avg1d, delt, dx, dxsq, rhoe, rhoed, th0, th0d, xlv, cp, rd&
&  , rv, g, ep2, svp1, svp2, svp3, svpt0, dqvdt, dqvdtd, dthdt, dthdtd, &
&  pprate, pprated, nca, ntst, cutop, cubot, ids, ide, jds, jde, kds, kde&
&  , ims, ime, jms, jme, kms, kme, its, ite, jts, jte, kts, kte)
  IMPLICIT NONE
!-----------------------------------------------------------------------
!-----------------------------------------------------------
  INTEGER, INTENT(IN) :: ids, ide, jds, jde, kds, kde, ims, ime, jms, &
&  jme, kms, kme, its, ite, jts, jte, kts, kte, i, j, ntst
!
  REAL, DIMENSION(kts:kte), INTENT(IN) :: u0, v0, t0, th0, qv0, p0, rhoe&
&  , dzq, z, w0avg1d
  REAL, DIMENSION(kts:kte), INTENT(IN) :: t0d, th0d, qv0d, p0d, rhoed, &
&  dzqd, zd
!
  REAL, INTENT(IN) :: delt, dx, dxsq
!
  REAL, INTENT(IN) :: xlv, cp, rd, rv, g
  REAL, INTENT(IN) :: ep2, svp1, svp2, svp3, svpt0
!
  REAL, DIMENSION(kts:kte), INTENT(INOUT) :: dqvdt, dthdt
  REAL, DIMENSION(kts:kte), INTENT(INOUT) :: dqvdtd, dthdtd
  REAL, DIMENSION(ims:ime, jms:jme), INTENT(INOUT) :: nca
  REAL, DIMENSION(ims:ime, jms:jme), INTENT(OUT) :: cubot, cutop
  REAL, INTENT(OUT) :: pprate
  REAL, INTENT(OUT) :: pprated
!
!...DEFINE LOCAL VARIABLES...
!
  REAL, DIMENSION(kts:kte) :: cond, h, hs, qs, x
  REAL, DIMENSION(kts:kte) :: condd, hd, hsd, qsd, xd
  REAL :: buoy, cape, cin, dh, dq, dt, dtm, ep, es, evap, hp, mp, qp, &
&  qsp, rrk, rrkp, tadp, tdp, zb, zg, zi, zt
  REAL :: dhd, dqd, dtd, esd, evapd, hpd, mpd, qpd, qspd, rrkd, rrkpd
  INTEGER :: ipos, isat, k, kb, ki, kt
  REAL :: arg1
  REAL :: arg1d
  REAL :: arg10
  REAL :: arg2
  REAL :: abs1
  REAL :: abs0
  hd = 0.0_8
  qsd = 0.0_8
  xd = 0.0_8
  hsd = 0.0_8
!
!...DEFINE PROFILES
  DO k=kts,kte
    hd(k) = cp*t0d(k) + g*zd(k) + xlv*qv0d(k)
    h(k) = cp*t0(k) + g*z(k) + xlv*qv0(k)
    arg1d = (svp2*t0d(k)*(t0(k)-svp3)-svp2*(t0(k)-svpt0)*t0d(k))/(t0(k)-&
&      svp3)**2
    arg1 = svp2*(t0(k)-svpt0)/(t0(k)-svp3)
    esd = 1000.*svp1*arg1d*EXP(arg1)
    es = 1000.*svp1*EXP(arg1)
    qsd(k) = (ep2*esd*(p0(k)-es)-ep2*es*(p0d(k)-esd))/(p0(k)-es)**2
    qs(k) = ep2*es/(p0(k)-es)
    hsd(k) = cp*t0d(k) + g*zd(k) + xlv*qsd(k)
    hs(k) = cp*t0(k) + g*z(k) + xlv*qs(k)
    xd(k) = (xlv**2*qsd(k)*cp*rv*t0(k)**2-xlv**2*qs(k)*cp*rv*(t0d(k)*t0(&
&      k)+t0(k)*t0d(k)))/(cp*rv*t0(k)*t0(k))**2
    x(k) = xlv*xlv*qs(k)/(cp*rv*t0(k)*t0(k))
    dthdtd(k) = 0.0_8
    dthdt(k) = 0.
    dqvdtd(k) = 0.0_8
    dqvdt(k) = 0.
  END DO
  pprate = 0.
  zg = z(1) - 0.5*dzq(1)
  pprated = 0.0_8
!
!...LOOP OVER PARCELS
loop_origin:DO ki=kts,kte
    hpd = hd(ki)
    hp = h(ki)
    qpd = qv0d(ki)
    qp = qv0(ki)
    mpd = alpha*(rhoed(ki)*dzq(ki)+rhoe(ki)*dzqd(ki))
    mp = alpha*rhoe(ki)*dzq(ki)
    zi = z(ki)
    buoy = 0.
    cape = 0.
    cin = 0.
    dtm = 0.
    isat = 0
    ipos = 0
    kt = 0
    kb = 0
    zt = 0.
    zb = 0.
    cond = 0.
    condd = 0.0_8
!
!...LIFT PARCEL
loop_lift:DO k=ki+1,kte
      tadp = t0(ki) + g/cp*(z(ki)-z(k))
      ep = p0(k)*qv0(ki)/(ep2+qv0(ki))
      arg10 = 0.001*ep/svp1
      arg2 = 0.001*ep/svp1
      tdp = (svpt0-svp3/svp2*ALOG(arg10))/(1.-1./svp2*ALOG(arg2))
      IF (tadp .GE. tdp) THEN
!         unsaturated
        IF (isat .EQ. 1) PRINT*, i, j, &
&                         'sounding warning: unsat above sat'
        dt = tadp - t0(k)
        condd(k) = 0.0_8
        cond(k) = 0.
      ELSE
!         saturated
        IF (isat .EQ. 0) THEN
          kb = k
          zb = z(k) - 0.5*dzq(k)
        END IF
        isat = 1
        dhd = hpd - hsd(k)
        dh = hp - hs(k)
        dt = dh/cp/(1.+x(k))
        qspd = qsd(k) + ((dhd*x(k)/xlv+dh*xd(k)/xlv)*(1.+x(k))-dh*x(k)*&
&          xd(k)/xlv)/(1.+x(k))**2
        qsp = qs(k) + dh/xlv*x(k)/(1.+x(k))
!...CONDENSATE PRODUCED
        condd(k) = mpd*(qp-qsp) + mp*(qpd-qspd)
        cond(k) = mp*(qp-qsp)
        qpd = qspd
        qp = qsp
      END IF
      buoy = buoy + g*dt*dzq(k)/t0(k)
      IF (cape .LT. buoy) THEN
        cape = buoy
      ELSE
        cape = cape
      END IF
      IF (buoy .GE. cincap) THEN
        IF (cin .GT. buoy) THEN
          cin = buoy
        ELSE
          cin = cin
        END IF
      END IF
      IF (dt .GE. 0.) THEN
        kt = k
        zt = z(k) + 0.5*dzq(k)
      ELSE IF (dt .LT. 0. .AND. dtm .GE. 0.) THEN
        IF (dt .GE. 0.) THEN
          abs0 = dt
        ELSE
          abs0 = -dt
        END IF
        IF (dtm .GE. 0.) THEN
          abs1 = dtm
        ELSE
          abs1 = -dtm
        END IF
! cloud top is level closest to parcel temperature
        IF (abs0 .LT. abs1) THEN
          kt = k
          zt = z(k) + 0.5*dzq(k)
        END IF
      END IF
      dtm = dt
! continue lifting until buoyancy is gone
      IF (buoy .LT. cincap) THEN
        GOTO 100
      ELSE
        IF (buoy .GT. 0.) ipos = 1
!         positive area detected
        IF (k .EQ. 1) THEN
          kt = k
          zt = z(k) + 0.5*dzq(k)
          zi = z(ki)
          PRINT*, 'sounding warning: cloud top at model top'
        END IF
      END IF
    END DO loop_lift
!
!...CHECK FOR CLOUD
 100 IF (isat .NE. 0) THEN
!       no cloud from lifting - no convection
      IF (zt - zb .GT. dpthmin) THEN
!       not more than one cloud level - no convection
        IF (ipos .NE. 0) THEN
!       no buoyancy in cloud - no convection
          IF (cape .GT. capemin) THEN
!       not enough cape
!
!...IF CHECK FOR CLOUD SUCCESSFUL
!
!...DETRAINMENT
            dhd = hpd - hsd(kt)
            dh = hp - hs(kt)
            dtd = (dhd*(1.+x(kt))/cp-dh*xd(kt)/cp)/(1.+x(kt))**2
            dt = dh/cp/(1.+x(kt))
            dqd = qsd(kt) + ((dhd*x(kt)/xlv+dh*xd(kt)/xlv)*(1.+x(kt))-dh&
&              *x(kt)*xd(kt)/xlv)/(1.+x(kt))**2 - qv0d(kt)
            dq = qs(kt) + dh/xlv*x(kt)/(1.+x(kt)) - qv0(kt)
            dthdtd(kt) = dthdtd(kt) + (((mpd*dt+mp*dtd)*th0(kt)/t0(kt)+&
&              mp*dt*(th0d(kt)*t0(kt)-th0(kt)*t0d(kt))/t0(kt)**2)*rhoe(kt&
&              )*dzq(kt)-mp*th0(kt)*dt*(rhoed(kt)*dzq(kt)+rhoe(kt)*dzqd(&
&              kt))/t0(kt))/(rhoe(kt)*dzq(kt))**2
            dthdt(kt) = dthdt(kt) + mp*(th0(kt)/t0(kt))*dt/(rhoe(kt)*dzq&
&              (kt))
            dqvdtd(kt) = dqvdtd(kt) + ((mpd*dq+mp*dqd)*rhoe(kt)*dzq(kt)-&
&              mp*dq*(rhoed(kt)*dzq(kt)+rhoe(kt)*dzqd(kt)))/(rhoe(kt)*dzq&
&              (kt))**2
            dqvdt(kt) = dqvdt(kt) + mp*dq/(rhoe(kt)*dzq(kt))
!
!...SUBSIDENCE
loop_subsidence:DO k=kt-1,ki,-1
              dthdtd(k) = dthdtd(k) + ((mpd*(th0(k+1)-th0(k))+mp*(th0d(k&
&                +1)-th0d(k)))*rhoe(k)*dzq(k)-mp*(th0(k+1)-th0(k))*(rhoed&
&                (k)*dzq(k)+rhoe(k)*dzqd(k)))/(rhoe(k)*dzq(k))**2
              dthdt(k) = dthdt(k) + mp*(th0(k+1)-th0(k))/(rhoe(k)*dzq(k)&
&                )
              dqvdtd(k) = dqvdtd(k) + ((mpd*(qv0(k+1)-qv0(k))+mp*(qv0d(k&
&                +1)-qv0d(k)))*rhoe(k)*dzq(k)-mp*(qv0(k+1)-qv0(k))*(rhoed&
&                (k)*dzq(k)+rhoe(k)*dzqd(k)))/(rhoe(k)*dzq(k))**2
              dqvdt(k) = dqvdt(k) + mp*(qv0(k+1)-qv0(k))/(rhoe(k)*dzq(k)&
&                )
            END DO loop_subsidence
!
!...RAINFALL AND EVAPORATION
            rrkp = 0.
            rrkpd = 0.0_8
loop_rainfall:DO k=kt,1,-1
              rrkd = rrkpd + condd(k)
              rrk = rrkp + cond(k)
              evapd = eps*((dzqd(k)*rrkp+dzq(k)*rrkpd)*(qs(k)-qv0(k))/&
&                vfall+dzq(k)*rrkp*(qsd(k)-qv0d(k))/vfall)
              evap = dzq(k)*rrkp/vfall*eps*(qs(k)-qv0(k))
! restrict evap to below cloud base
              IF (k .GE. kb) THEN
                evap = 0.
                evapd = 0.0_8
              END IF
              IF (rrk .GT. evap) THEN
                evap = evap
              ELSE
                evapd = rrkd
                evap = rrk
              END IF
              rrkd = rrkd - evapd
              rrk = rrk - evap
              dqvdtd(k) = dqvdtd(k) + (evapd*rhoe(k)*dzq(k)-evap*(rhoed(&
&                k)*dzq(k)+rhoe(k)*dzqd(k)))/(rhoe(k)*dzq(k))**2
              dqvdt(k) = dqvdt(k) + evap/(rhoe(k)*dzq(k))
              dthdtd(k) = dthdtd(k) - (xlv*((th0d(kt)*t0(kt)-th0(kt)*t0d&
&                (kt))*evap/t0(kt)**2+th0(kt)*evapd/t0(kt))*rhoe(k)*dzq(k&
&                )/cp-xlv*th0(kt)*evap*(rhoed(k)*dzq(k)+rhoe(k)*dzqd(k))/&
&                (cp*t0(kt)))/(rhoe(k)*dzq(k))**2
              dthdt(k) = dthdt(k) - xlv/cp*(th0(kt)/t0(kt))*evap/(rhoe(k&
&                )*dzq(k))
              rrkpd = rrkd
              rrkp = rrk
            END DO loop_rainfall
            pprated = pprated + rrkpd
            pprate = pprate + rrkp
          END IF
        END IF
      END IF
    END IF
  END DO loop_origin
END SUBROUTINE DUCU1D_D

!        Generated by TAPENADE     (INRIA, Tropics team)
!  Tapenade 3.6 (r4343) - 10 Feb 2012 10:52
!
!  Differentiation of ducuinit in forward (tangent) mode:
!   variations   of useful results: rqccuten rthcuten w0avg rqvcuten
!   with respect to varying inputs: rqccuten rthcuten w0avg rqvcuten
SUBROUTINE DUCUINIT_D(rthcuten, rthcutend, rqvcuten, rqvcutend, rqccuten&
&  , rqccutend, rqrcuten, rqicuten, rqscuten, nca, w0avg, w0avgd, p_qc, &
&  p_qr, svp1, svp2, svp3, svpt0, p_first_scalar, restart, &
&  allowed_to_read, ids, ide, jds, jde, kds, kde, ims, ime, jms, jme, kms&
&  , kme, its, ite, jts, jte, kts, kte)
  IMPLICIT NONE
!--------------------------------------------------------------------
  LOGICAL, INTENT(IN) :: restart, allowed_to_read
  INTEGER, INTENT(IN) :: ids, ide, jds, jde, kds, kde, ims, ime, jms, &
&  jme, kms, kme, its, ite, jts, jte, kts, kte
  INTEGER, INTENT(IN) :: p_qc, p_qr, p_first_scalar
  REAL, DIMENSION(ims:ime, kms:kme, jms:jme), INTENT(OUT) :: rthcuten, &
&  rqvcuten, rqccuten, rqrcuten, rqicuten, rqscuten
  REAL, DIMENSION(ims:ime, kms:kme, jms:jme), INTENT(OUT) :: rthcutend, &
&  rqvcutend, rqccutend
  REAL, DIMENSION(ims:ime, kms:kme, jms:jme), INTENT(OUT) :: w0avg
  REAL, DIMENSION(ims:ime, kms:kme, jms:jme), INTENT(OUT) :: w0avgd
  REAL, DIMENSION(ims:ime, jms:jme), INTENT(INOUT) :: nca
  INTEGER :: i, j, k, itf, jtf, ktf
  REAL, INTENT(IN) :: svp1, svp2, svp3, svpt0
  IF (jte .GT. jde - 1) THEN
    jtf = jde - 1
  ELSE
    jtf = jte
  END IF
  IF (kte .GT. kde - 1) THEN
    ktf = kde - 1
  ELSE
    ktf = kte
  END IF
  IF (ite .GT. ide - 1) THEN
    itf = ide - 1
  ELSE
    itf = ite
  END IF
  IF (.NOT.restart) THEN
    DO j=jts,jtf
      DO k=kts,ktf
        DO i=its,itf
          rthcutend(i, k, j) = 0.0
          rthcuten(i, k, j) = 0.
          rqvcutend(i, k, j) = 0.0
          rqvcuten(i, k, j) = 0.
        END DO
      END DO
    END DO
    IF (p_qc .GE. p_first_scalar) THEN
      DO j=jts,jtf
        DO k=kts,ktf
          DO i=its,itf
            rqccutend(i, k, j) = 0.0
            rqccuten(i, k, j) = 0.
          END DO
        END DO
      END DO
    END IF
    DO j=jts,jtf
      DO i=its,itf
        nca(i, j) = -100.
      END DO
    END DO
    DO j=jts,jtf
      DO k=kts,ktf
        DO i=its,itf
          w0avgd(i, k, j) = 0.0
          w0avg(i, k, j) = 0.
        END DO
      END DO
    END DO
  END IF
END SUBROUTINE DUCUINIT_D

END MODULE g_module_cu_du