cloud_tl.F90

Go to the documentation of this file.

MODULE cloud_tl

USE cloud, ONLY: pdf_width

IMPLICIT NONE

PRIVATE
PUBLIC cloud_driver_d, ls_cloud_d

CONTAINS
!        Generated by TAPENADE     (INRIA, Tropics team)
!  Tapenade 3.9 (r5096) - 24 Feb 2014 16:53
!
!  Differentiation of cloud_driver in forward (tangent) mode:
!   variations   of useful results: th qi_con q cf_ls cf_con ql_ls
!                ql_con qi_ls
!   with respect to varying inputs: th qi_con cnv_prc3 q cnv_dqldt
!                cf_ls cnv_updf cf_con cnv_mfd ql_ls ql_con qi_ls
!   RW status of diff variables: th:in-out qi_con:in-out cnv_prc3:in
!                q:in-out cnv_dqldt:in cf_ls:in-out cnv_updf:in
!                cf_con:in-out cnv_mfd:in ql_ls:in-out ql_con:in-out
!                qi_ls:in-out
SUBROUTINE cloud_driver_d(dt, im, jm, lm, th, thd, q, qd, ple, cnv_dqldt&
& , cnv_dqldtd, cnv_mfd, cnv_mfdd, cnv_prc3, cnv_prc3d, cnv_updf, &
& cnv_updfd, qi_ls, qi_lsd, ql_ls, ql_lsd, qi_con, qi_cond, ql_con, &
& ql_cond, cf_ls, cf_lsd, cf_con, cf_cond, frland, physparams, estblx, &
& khu, khl, cons_runiv, cons_kappa, cons_airmw, cons_h2omw, cons_grav, &
& cons_alhl, cons_alhf, cons_pi, cons_rgas, cons_cp, cons_vireps, &
& cons_alhs, cons_tice, cons_rvap, cons_p00, do_moist_physics)
  IMPLICIT NONE
!INPUTS
  INTEGER, INTENT(IN) :: im, jm, lm, do_moist_physics
  REAL*8, INTENT(IN) :: dt, frland(im, jm), physparams(:)
  REAL*8, DIMENSION(im, jm, lm), INTENT(IN) :: cnv_dqldt, cnv_mfd, &
& cnv_updf, cnv_prc3
  REAL*8, DIMENSION(im, jm, lm), INTENT(IN) :: cnv_dqldtd, cnv_mfdd, &
& cnv_updfd, cnv_prc3d
  REAL*8, INTENT(IN) :: estblx(:)
  REAL*8, DIMENSION(im, jm, 0:lm), INTENT(IN) :: ple
  INTEGER, DIMENSION(im, jm), INTENT(IN) :: khu, khl
!MAPL_CONSTANTS REDEFINED FOR USE IN AUTODIFF TOOL
  REAL*8, INTENT(IN) :: cons_runiv, cons_kappa, cons_airmw
  REAL*8, INTENT(IN) :: cons_h2omw, cons_grav, cons_alhl
  REAL*8, INTENT(IN) :: cons_alhf, cons_pi, cons_rgas
  REAL*8, INTENT(IN) :: cons_cp, cons_vireps, cons_alhs
  REAL*8, INTENT(IN) :: cons_tice, cons_rvap, cons_p00
!PROGNOSTICS
  REAL*8, DIMENSION(im, jm, lm), INTENT(INOUT) :: th, q
  REAL*8, DIMENSION(im, jm, lm), INTENT(INOUT) :: thd, qd
  REAL*8, DIMENSION(im, jm, lm), INTENT(INOUT) :: qi_ls, ql_ls, qi_con, &
& ql_con, cf_con, cf_ls
  REAL*8, DIMENSION(im, jm, lm), INTENT(INOUT) :: qi_lsd, ql_lsd, &
& qi_cond, ql_cond, cf_cond, cf_lsd
!OUTPUTS (DIAGNOSTICS)
!LOCALS
  INTEGER :: i, j, k, l, ktop
  REAL*8, DIMENSION(im, jm, lm) :: ph, pih, mass, imass, t, dzet, qddf3&
& , rh, dp, dm
  REAL*8, DIMENSION(im, jm, lm) :: td, dzetd, qddf3d
  REAL*8, DIMENSION(im, jm, lm + 1) :: zet
  REAL*8, DIMENSION(im, jm, lm+1) :: zetd
  REAL*8, DIMENSION(im, jm, 0:lm) :: p, pi
  REAL*8, DIMENSION(im, jm, lm) :: qs, dqsdt, dqs
  REAL*8, DIMENSION(im, jm, lm) :: qsd, dqsdtd, dqsd
  REAL*8, DIMENSION(im, jm) :: vmip
  REAL*8, DIMENSION(im, jm) :: vmipd
!Precip amounts and fall rate
  REAL*8 :: cf_tot
  REAL*8 :: cf_totd
  REAL*8 :: alpha, alhx3, rhcrit
  REAL*8 :: alhx3d
!Microphyiscal constants
  REAL*8 :: aa, bb
  REAL*8 :: aad, bbd
  REAL*8, PARAMETER :: pi_0=4.*atan(1.)
!Density of liquid water in kg/m^3
  REAL*8, PARAMETER :: rho_w=1.0e3
  REAL*8 :: t_ice_max
!PHYSPARAMS constants
  REAL*8 :: cnv_beta, anv_beta, ls_beta, rh00, c_00, lwcrit, c_acc, &
& c_ev_r, c_ev_s, cldvol2frc
  REAL*8 :: rhsup_ice, shr_evap_fac, min_cld_water, cld_evp_eff, &
& ls_sdqv2, ls_sdqv3, ls_sdqvt1
  REAL*8 :: anv_sdqv2, anv_sdqv3, anv_sdqvt1, anv_to_ls, n_warm, n_ice, &
& n_anvil, n_pbl
  REAL*8 :: anv_icefall_c, ls_icefall_c, revap_off_p, cnvenvfc, wrhodep&
& , t_ice_all, cnviceparam
  REAL*8 :: cnvddrfc, anvddrfc, lsddrfc, minrhcrit, maxrhcrit, &
& turnrhcrit, maxrhcritland
  REAL*8 :: min_rl, min_ri, max_rl, max_ri, ri_anv
  INTEGER :: nsmax, disable_rad, icefrpwr, tanhrhcrit, fr_ls_wat, &
& fr_ls_ice, fr_an_wat, fr_an_ice, pdfflag
  REAL*8 :: lsenvfc, anvenvfc
  REAL*8 :: qrn_cu, qsn_cu, qrn_an, qsn_an, qrn_ls, qsn_ls, qrn_cu_1d
  REAL*8 :: qsn_cud, qrn_and, qsn_and, qrn_lsd, qsn_lsd, qrn_cu_1dd
  REAL*8 :: qt_tmpi_1, qt_tmpi_2, qlt_tmp, qit_tmp
  REAL*8 :: qt_tmpi_1d, qt_tmpi_2d, qlt_tmpd, qit_tmpd
  REAL*8 :: prn_above_cu_new, prn_above_an_new, prn_above_ls_new
  REAL*8 :: prn_above_cu_newd, prn_above_an_newd, prn_above_ls_newd
  REAL*8 :: prn_above_cu_old, prn_above_an_old, prn_above_ls_old
  REAL*8 :: prn_above_cu_oldd, prn_above_an_oldd, prn_above_ls_oldd
  REAL*8 :: psn_above_cu_new, psn_above_an_new, psn_above_ls_new
  REAL*8 :: psn_above_cu_newd, psn_above_an_newd, psn_above_ls_newd
  REAL*8 :: psn_above_cu_old, psn_above_an_old, psn_above_ls_old
  REAL*8 :: psn_above_cu_oldd, psn_above_an_oldd, psn_above_ls_oldd
  REAL*8 :: evap_dd_cu_above_new, evap_dd_an_above_new, &
& evap_dd_ls_above_new
  REAL*8 :: evap_dd_cu_above_newd, evap_dd_an_above_newd, &
& evap_dd_ls_above_newd
  REAL*8 :: evap_dd_cu_above_old, evap_dd_an_above_old, &
& evap_dd_ls_above_old
  REAL*8 :: evap_dd_cu_above_oldd, evap_dd_an_above_oldd, &
& evap_dd_ls_above_oldd
  REAL*8 :: subl_dd_cu_above_new, subl_dd_an_above_new, &
& subl_dd_ls_above_new
  REAL*8 :: subl_dd_cu_above_newd, subl_dd_an_above_newd, &
& subl_dd_ls_above_newd
  REAL*8 :: subl_dd_cu_above_old, subl_dd_an_above_old, &
& subl_dd_ls_above_old
  REAL*8 :: subl_dd_cu_above_oldd, subl_dd_an_above_oldd, &
& subl_dd_ls_above_oldd
  REAL*8 :: area_ls_prc1, area_upd_prc1, area_anv_prc1
  REAL*8 :: area_ls_prc1d, area_upd_prc1d, area_anv_prc1d
  REAL*8 :: tot_prec_upd, tot_prec_anv, tot_prec_ls, area_upd_prc, &
& area_anv_prc, area_ls_prc
  REAL*8 :: tot_prec_updd, tot_prec_anvd, tot_prec_lsd, area_upd_prcd, &
& area_anv_prcd, area_ls_prcd
  REAL*8 :: qtmp2
  REAL*8 :: rhexcess, tpw, negtpw
  REAL*8 :: tpwd, negtpwd
  INTEGER :: cloud_pertmod
  INTRINSIC atan
  INTRINSIC int
  INTRINSIC sum
  INTRINSIC max
  REAL*8, DIMENSION(im, jm, 0:lm) :: pwx1
  REAL*8 :: pwy1
  REAL*8, DIMENSION(im, jm, lm) :: pwx10

  !LS_CLOUD FILTERING
  integer :: ii
  real(8) :: xx(8) 
  real(8) :: ttraj, qtraj, qi_lstraj, qi_contraj, ql_lstraj, ql_contraj, cf_lstraj, cf_contraj, phtraj 
  real(8) :: tpert, qpert, qi_lspert, qi_conpert, ql_lspert, ql_conpert, cf_lspert, cf_conpert
  real(8) :: jacobian(8,8), a(8,8)

  !Eigenvalue computation
  integer, parameter :: n1  = 8
  integer, parameter :: lda = 8
  REAL(8) :: wr(n1), wi(n1)
  INTEGER, PARAMETER :: ldvl = n1 !Left vecotrs
  REAL(8) :: vl(ldvl,n1)
  INTEGER, PARAMETER :: ldvr = n1 !Right vectors
  REAL(8) :: vr(ldvr,n1)
  INTEGER :: info, lwork
  INTEGER, PARAMETER :: lwmax = 10000
  REAL(8) :: work(lwmax)
  EXTERNAL :: dgeev
  REAL*8 :: maxeval

  !Total Filtering
  real(8) :: totfilt_t, totfilt_ql, totfilt_qi
  real(8) :: t_p_preall, ql_ls_p_preall, ql_con_p_preall, qi_ls_p_preall, qi_con_p_preall

  !Sink Filtering
  real(8) :: sinkfilt_ql, sinkfilt_qi, sinkfilt_cf
  real(8) :: t_p_presink, q_p_presink
  real(8) ::  ql_ls_p_presink, ql_con_p_presink
  real(8) ::  qi_ls_p_presink, qi_con_p_presink
  !real(8) ::  cf_con_p_presink

!Highest level of calculations
  ktop = 30
!Get Constants from CLOUDPARAMS
! Area factor for convective rain showers (non-dim)
  cnv_beta = physparams(1)
! Area factor for anvil rain showers (non-dim)
  anv_beta = physparams(2)
! Area factor for Large Scale rain showers (non-dim)
  ls_beta = physparams(3)
! Critical relative humidity
  rh00 = physparams(4)
  c_00 = physparams(5)
  lwcrit = physparams(6)
  c_acc = physparams(7)
  c_ev_r = physparams(8)
  c_ev_s = physparams(56)
  cldvol2frc = physparams(9)
  rhsup_ice = physparams(10)
  shr_evap_fac = physparams(11)
  min_cld_water = physparams(12)
  cld_evp_eff = physparams(13)
  nsmax = int(physparams(14))
  ls_sdqv2 = physparams(15)
  ls_sdqv3 = physparams(16)
  ls_sdqvt1 = physparams(17)
  anv_sdqv2 = physparams(18)
  anv_sdqv3 = physparams(19)
  anv_sdqvt1 = physparams(20)
  anv_to_ls = physparams(21)
  n_warm = physparams(22)
  n_ice = physparams(23)
  n_anvil = physparams(24)
  n_pbl = physparams(25)
  disable_rad = int(physparams(26))
  anv_icefall_c = physparams(28)
  ls_icefall_c = physparams(29)
  revap_off_p = physparams(30)
  cnvenvfc = physparams(31)
  wrhodep = physparams(32)
  t_ice_all = physparams(33) + cons_tice
  cnviceparam = physparams(34)
  icefrpwr = int(physparams(35) + .001)
  cnvddrfc = physparams(36)
  anvddrfc = physparams(37)
  lsddrfc = physparams(38)
  tanhrhcrit = int(physparams(41))
  minrhcrit = physparams(42)
  maxrhcrit = physparams(43)
  turnrhcrit = physparams(45)
  maxrhcritland = physparams(46)
  fr_ls_wat = int(physparams(47))
  fr_ls_ice = int(physparams(48))
  fr_an_wat = int(physparams(49))
  fr_an_ice = int(physparams(50))
  min_rl = physparams(51)
  min_ri = physparams(52)
  max_rl = physparams(53)
  max_ri = physparams(54)
  ri_anv = physparams(55)
  pdfflag = int(physparams(57))
  t_ice_max = cons_tice
!Initialize the saving of downdraft values.
  prn_above_cu_new = 0.
  prn_above_an_new = 0.
  prn_above_ls_new = 0.
  psn_above_cu_new = 0.
  psn_above_an_new = 0.
  psn_above_ls_new = 0.
  evap_dd_cu_above_new = 0.
  evap_dd_an_above_new = 0.
  evap_dd_ls_above_new = 0.
  subl_dd_cu_above_new = 0.
  subl_dd_an_above_new = 0.
  subl_dd_ls_above_new = 0.
!Convert to hPa and average pressure to temperature levels
  p = ple*0.01
  ph = 0.5*(p(:, :, 0:lm-1)+p(:, :, 1:lm))
!Calculate Exner Pressure at temperature levels
  pwx1 = p/1000.
  pwy1 = cons_rgas/cons_cp
  pi = pwx1**pwy1
  pwx10 = ph/1000.
  pwy1 = cons_rgas/cons_cp
  pih = pwx10**pwy1
!Calculate temperature
  td = pih*thd
  t = th*pih
!Compute QS and DQSDT
  qsd = 0.0_8
  dqsdtd = 0.0_8
  CALL dqsat_bac_d(dqsdt, dqsdtd, qs, qsd, t, td, ph, im, jm, lm, estblx&
&            , cons_h2omw, cons_airmw)
!Relative humidity
  rh = q/qs
!Compute layer mass and 1/mass
  mass = (p(:, :, 1:lm)-p(:, :, 0:lm-1))*100./cons_grav
  imass = 1/mass
!Level thickness
  dzetd(:, :, 1:lm) = (pi(:, :, 1:lm)-pi(:, :, 0:lm-1))*cons_cp*thd(:, :&
&   , 1:lm)/cons_grav
  dzet(:, :, 1:lm) = th(:, :, 1:lm)*(pi(:, :, 1:lm)-pi(:, :, 0:lm-1))*&
&   cons_cp/cons_grav
!Level heights
  zetd(:, :, lm+1) = 0.0_8
  zet(:, :, lm+1) = 0.0
  zetd = 0.0_8
  DO k=lm,1,-1
    zetd(:, :, k) = zetd(:, :, k+1) + dzetd(:, :, k)
    zet(:, :, k) = zet(:, :, k+1) + dzet(:, :, k)
  END DO
  qddf3d = 0.0_8
  WHERE (zet(:, :, 1:lm) .LT. 3000.) 
    qddf3d = -(mass*(zetd(:, :, 1:lm)*zet(:, :, 1:lm)+(zet(:, :, 1:lm)-&
&     3000.)*zetd(:, :, 1:lm)))
    qddf3 = -((zet(:, :, 1:lm)-3000.)*zet(:, :, 1:lm)*mass)
  ELSEWHERE
    qddf3d = 0.0_8
    qddf3 = 0.
  END WHERE
  vmipd = 0.0_8
  DO i=1,im
    DO j=1,jm
      vmipd(i, j) = sum(qddf3d(i, j, :))
      vmip(i, j) = sum(qddf3(i, j, :))
    END DO
  END DO
  DO k=1,lm
    qddf3d(:, :, k) = (qddf3d(:, :, k)*vmip-qddf3(:, :, k)*vmipd)/vmip**&
&     2
    qddf3(:, :, k) = qddf3(:, :, k)/vmip
  END DO
!Pressure and mass thickness for use in cleanup.
  dp = ple(:, :, 1:lm) - ple(:, :, 0:lm-1)
  dm = dp*(1./cons_grav)
  prn_above_cu_newd = 0.0_8
  psn_above_cu_newd = 0.0_8
  evap_dd_an_above_newd = 0.0_8
  area_upd_prcd = 0.0_8
  prn_above_ls_newd = 0.0_8
  evap_dd_ls_above_newd = 0.0_8
  evap_dd_cu_above_newd = 0.0_8
  psn_above_ls_newd = 0.0_8
  tot_prec_lsd = 0.0_8
  area_ls_prcd = 0.0_8
  tot_prec_updd = 0.0_8
  prn_above_an_newd = 0.0_8
  area_anv_prcd = 0.0_8
  psn_above_an_newd = 0.0_8
  alhx3d = 0.0_8
  subl_dd_an_above_newd = 0.0_8
  tot_prec_anvd = 0.0_8
  subl_dd_ls_above_newd = 0.0_8
  subl_dd_cu_above_newd = 0.0_8
!Begin loop over all grid boxes.
  DO i=1,im
    DO j=1,jm
      DO k=ktop,lm

        !Save the inputs to the scheme for filtering
        t_p_preall = td(i, j, k)
        ql_ls_p_preall = ql_lsd(i, j, k)
        ql_con_p_preall = ql_cond(i, j, k)
        qi_ls_p_preall = qi_lsd(i, j, k)
        qi_con_p_preall = qi_cond(i, j, k)

        IF (k .EQ. ktop) THEN
          tot_prec_upd = 0.
          tot_prec_anv = 0.
          tot_prec_ls = 0.
          area_upd_prc = 0.
          area_anv_prc = 0.
          area_ls_prc = 0.
          area_upd_prcd = 0.0_8
          tot_prec_lsd = 0.0_8
          area_ls_prcd = 0.0_8
          tot_prec_updd = 0.0_8
          area_anv_prcd = 0.0_8
          tot_prec_anvd = 0.0_8
        END IF
!Initialize precips, except QRN_CU which comes from RAS 
        qrn_ls = 0.
        qrn_an = 0.
        qrn_cu_1d = 0.
        qsn_ls = 0.
        qsn_an = 0.
        qsn_cu = 0.
!Ras Rain         
        qrn_cu_1dd = cnv_prc3d(i, j, k)
        qrn_cu_1d = cnv_prc3(i, j, k)
!Tidy up where fractions or cloud is too low
        CALL cloud_tidy_d(q(i, j, k), qd(i, j, k), t(i, j, k), td(i, j, &
&                   k), ql_ls(i, j, k), ql_lsd(i, j, k), qi_ls(i, j, k)&
&                   , qi_lsd(i, j, k), cf_ls(i, j, k), cf_lsd(i, j, k), &
&                   ql_con(i, j, k), ql_cond(i, j, k), qi_con(i, j, k), &
&                   qi_cond(i, j, k), cf_con(i, j, k), cf_cond(i, j, k)&
&                   , cons_alhl, cons_alhs, cons_cp)
!Phase changes for large scale cloud.
        CALL meltfreeze_d(dt, t(i, j, k), td(i, j, k), ql_ls(i, j, k), &
&                   ql_lsd(i, j, k), qi_ls(i, j, k), qi_lsd(i, j, k), &
&                   t_ice_all, t_ice_max, icefrpwr, cons_alhl, cons_alhs&
&                   , cons_cp)
!Phase changes for convective cloud.
        CALL meltfreeze_d(dt, t(i, j, k), td(i, j, k), ql_con(i, j, k), &
&                   ql_cond(i, j, k), qi_con(i, j, k), qi_cond(i, j, k)&
&                   , t_ice_all, t_ice_max, icefrpwr, cons_alhl, &
&                   cons_alhs, cons_cp)
!STAGE 1 - Compute convective clouds from RAS diagnostics
        CALL convec_src_d(dt, mass(i, j, k), imass(i, j, k), t(i, j, k)&
&                   , td(i, j, k), q(i, j, k), qd(i, j, k), cnv_dqldt(i&
&                   , j, k), cnv_dqldtd(i, j, k), cnv_mfd(i, j, k), &
&                   cnv_mfdd(i, j, k), ql_con(i, j, k), ql_cond(i, j, k)&
&                   , qi_con(i, j, k), qi_cond(i, j, k), cf_con(i, j, k)&
&                   , cf_cond(i, j, k), qs(i, j, k), qsd(i, j, k), &
&                   cons_alhs, cons_alhl, cons_cp, t_ice_all, t_ice_max&
&                   , icefrpwr)
!STAGE 2a - Get PDF attributes
        CALL pdf_width(ph(i, j, k), frland(i, j), maxrhcrit, &
&                maxrhcritland, turnrhcrit, minrhcrit, pi_0, alpha)
        IF (alpha .LT. 1.0 - rh00) THEN
          alpha = 1.0 - rh00
        ELSE
          alpha = alpha
        END IF
        rhcrit = 1.0 - alpha

!STAGE 2b - Use PDF to compute large scale cloud effects and diagnostics,
!           also update convection clouds

        if (do_moist_physics == 1) then

           !For data assimilation purposes always employ the perturbation model for cloud fraction
           cloud_pertmod = 1

        elseif (do_moist_physics == 2) then

           !For longer forecasts only employ perturabtion model if filtering deems necessary
           cloud_pertmod = 0

           !COMPUTE THE JACOBIAN OF LS_CLOUD AND USE RESULTS TO FILTER 'BAD' POINTS 
           jacobian = 0.0
           DO ii = 1,8
              xx = 0.0
              xx(ii) = 1.0

              !Store seperate trajectory so as not to overwrite
              ttraj      = t(i,j,k)
              qtraj      = q(i,j,k)
              qi_lstraj  = qi_ls(i,j,k)
              qi_contraj = qi_con(i,j,k)
              ql_lstraj  = ql_ls(i,j,k)
              ql_contraj = ql_con(i,j,k)
              cf_lstraj  = cf_ls(i,j,k)
              cf_contraj = cf_con(i,j,k)
              phtraj     = ph(i,j,k)

              tpert      = xx(1)
              qpert      = xx(2)
              qi_lspert  = xx(3)
              qi_conpert = xx(4)
              ql_lspert  = xx(5)
              ql_conpert = xx(6)
              cf_lspert  = xx(7)
              cf_conpert = xx(8)

              CALL ls_cloud_d(dt, alpha, pdfflag, phtraj, ttraj, tpert, qtraj, qpert, &
                              ql_lstraj, ql_lspert, ql_contraj, ql_conpert, &
                              qi_lstraj, qi_lspert, qi_contraj, qi_conpert, &
                              cf_lstraj, cf_lspert, cf_contraj, cf_conpert, &
                              cons_alhl, cons_alhf, cons_alhs, cons_cp, cons_h2omw, cons_airmw, &
                              t_ice_all, t_ice_max, icefrpwr, estblx, 0, do_moist_physics)

              jacobian(1,ii) = tpert
              jacobian(2,ii) = qpert
              jacobian(3,ii) = qi_lspert
              jacobian(4,ii) = qi_conpert
              jacobian(5,ii) = ql_lspert
              jacobian(6,ii) = ql_conpert
              jacobian(7,ii) = cf_lspert
              jacobian(8,ii) = cf_conpert

           endDO

           !Compute eigenvalues of the Jacobian
           a = jacobian

           lwork = -1
           work = 0.0
           info = 0

           CALL dgeev( 'Vectors', 'Vectors', n1, a, lda, wr, wi, vl, ldvl, vr, ldvr, work, lwork, info )

           lwork = min( lwmax, int( work( 1 ) ) )
           CALL dgeev( 'Vectors', 'Vectors', n1, a, lda, wr, wi, vl, ldvl, vr, ldvr, work, lwork, info )

           maxeval = maxval(abs(wr))

           !Filter based on the eigenvalues
           if (maxeval > 1.001) then
              cloud_pertmod = 1
           endif

           !Filter based on the Jacobian values
           if ( ( jacobian(1,1) <  0.6     )   .or. &
                ( jacobian(2,1) >  0.75e-4 )   .or. &
                ( jacobian(5,1) < -0.75e-4 )   .or. &
                ( jacobian(7,1) < -1.10    ) ) then
              cloud_pertmod = 1
           endif

        endif

        CALL ls_cloud_d(dt, alpha, pdfflag, ph(i, j, k), t(i, j, k), td(&
&                 i, j, k), q(i, j, k), qd(i, j, k), ql_ls(i, j, k), &
&                 ql_lsd(i, j, k), ql_con(i, j, k), ql_cond(i, j, k), &
&                 qi_ls(i, j, k), qi_lsd(i, j, k), qi_con(i, j, k), &
&                 qi_cond(i, j, k), cf_ls(i, j, k), cf_lsd(i, j, k), &
&                 cf_con(i, j, k), cf_cond(i, j, k), cons_alhl, &
&                 cons_alhf, cons_alhs, cons_cp, cons_h2omw, cons_airmw&
&                 , t_ice_all, t_ice_max, icefrpwr, estblx, &
&                 cloud_pertmod, do_moist_physics)

        !SAVE PRESINKS INPUTS
        t_p_presink = td(i,j,k)
        q_p_presink = qd(i,j,k)
        qi_ls_p_presink = qi_lsd(i,j,k)
        qi_con_p_presink = qi_cond(i,j,k)
        ql_ls_p_presink = ql_lsd(i,j,k)
        ql_con_p_presink = ql_cond(i,j,k)
        !cf_con_p_presink = cf_cond(i,j,k)

!Clean up where too much overall cloud.
        cf_totd = cf_lsd(i, j, k) + cf_cond(i, j, k)
        cf_tot = cf_ls(i, j, k) + cf_con(i, j, k)
        IF (cf_tot .GT. 1.00) THEN
          cf_lsd(i, j, k) = cf_lsd(i, j, k)/cf_tot - cf_ls(i, j, k)*&
&           cf_totd/cf_tot**2
          cf_ls(i, j, k) = cf_ls(i, j, k)*(1.00/cf_tot)
          cf_cond(i, j, k) = cf_cond(i, j, k)/cf_tot - cf_con(i, j, k)*&
&           cf_totd/cf_tot**2
          cf_con(i, j, k) = cf_con(i, j, k)*(1.00/cf_tot)
        END IF
        cf_tot = cf_ls(i, j, k) + cf_con(i, j, k)
!STAGE 3 - Evap, Sublimation and Autoconversion
!Evaporation and sublimation of anvil cloud
        CALL evap_cnv_d(dt, rhcrit, ph(i, j, k), t(i, j, k), td(i, j, k)&
&                 , q(i, j, k), qd(i, j, k), ql_con(i, j, k), ql_cond(i&
&                 , j, k), qi_con(i, j, k), qi_cond(i, j, k), cf_con(i, &
&                 j, k), cf_cond(i, j, k), cf_ls(i, j, k), qs(i, j, k), &
&                 qsd(i, j, k), rho_w, cld_evp_eff, cons_h2omw, &
&                 cons_airmw, cons_alhl, cons_rvap, cons_rgas, cons_pi, &
&                 cons_cp)
        CALL subl_cnv_d(dt, rhcrit, ph(i, j, k), t(i, j, k), td(i, j, k)&
&                 , q(i, j, k), qd(i, j, k), ql_con(i, j, k), ql_cond(i&
&                 , j, k), qi_con(i, j, k), qi_cond(i, j, k), cf_con(i, &
&                 j, k), cf_cond(i, j, k), cf_ls(i, j, k), qs(i, j, k), &
&                 qsd(i, j, k), rho_w, cld_evp_eff, cons_h2omw, &
&                 cons_airmw, cons_alhl, cons_rvap, cons_rgas, cons_pi, &
&                 cons_cp, cons_alhs)
!Autoconversion
        CALL autoconversion_ls_d(dt, ql_ls(i, j, k), ql_lsd(i, j, k), &
&                          qrn_ls, qrn_lsd, t(i, j, k), td(i, j, k), ph(&
&                          i, j, k), cf_ls(i, j, k), cf_lsd(i, j, k), &
&                          ls_sdqv2, ls_sdqv3, ls_sdqvt1, c_00, lwcrit, &
&                          dzet(i, j, k))
        CALL autoconversion_cnv_d(dt, ql_con(i, j, k), ql_cond(i, j, k)&
&                           , qrn_an, qrn_and, t(i, j, k), td(i, j, k), &
&                           ph(i, j, k), cf_con(i, j, k), cf_cond(i, j, &
&                           k), anv_sdqv2, anv_sdqv3, anv_sdqvt1, c_00, &
&                           lwcrit, dzet(i, j, k))
!STAGE 4 - Fall and Re-evaporation of precip
        CALL ice_settlefall_cnv_d(wrhodep, qi_con(i, j, k), qi_cond(i, j&
&                           , k), ph(i, j, k), t(i, j, k), td(i, j, k), &
&                           cf_con(i, j, k), cf_cond(i, j, k), cons_rgas&
&                           , khu(i, j), khl(i, j), k, dt, dzet(i, j, k)&
&                           , dzetd(i, j, k), qsn_an, qsn_and, &
&                           anv_icefall_c)
        CALL ice_settlefall_ls_d(wrhodep, qi_ls(i, j, k), qi_lsd(i, j, k&
&                          ), ph(i, j, k), t(i, j, k), td(i, j, k), &
&                          cf_ls(i, j, k), cf_lsd(i, j, k), cons_rgas, &
&                          khu(i, j), khl(i, j), k, dt, dzet(i, j, k), &
&                          dzetd(i, j, k), qsn_ls, qsn_lsd, ls_icefall_c&
&                         )
!"Freeze" out any conv. precip, not done in RAS. This is
! precip w/ large particles, so freezing is strict.
        qtmp2 = 0.
        IF (t(i, j, k) .LT. cons_tice) THEN
          qtmp2 = qrn_cu_1d
          qsn_cud = qrn_cu_1dd
          qsn_cu = qrn_cu_1d
          qrn_cu_1d = 0.
          td(i, j, k) = td(i, j, k) + (cons_alhs-cons_alhl)*qsn_cud/&
&           cons_cp
          t(i, j, k) = t(i, j, k) + qsn_cu*(cons_alhs-cons_alhl)/cons_cp
          qrn_cu_1dd = 0.0_8
        ELSE
          qsn_cud = 0.0_8
        END IF
!Area
        area_ls_prc1 = 0.0
        area_upd_prc1 = 0.0
        area_anv_prc1 = 0.0
        tot_prec_updd = tot_prec_updd + mass(i, j, k)*(qrn_cu_1dd+&
&         qsn_cud)
        tot_prec_upd = tot_prec_upd + (qrn_cu_1d+qsn_cu)*mass(i, j, k)
        area_upd_prcd = area_upd_prcd + mass(i, j, k)*(cnv_updfd(i, j, k&
&         )*(qrn_cu_1d+qsn_cu)+cnv_updf(i, j, k)*(qrn_cu_1dd+qsn_cud))
        area_upd_prc = area_upd_prc + cnv_updf(i, j, k)*(qrn_cu_1d+&
&         qsn_cu)*mass(i, j, k)
        tot_prec_anvd = tot_prec_anvd + mass(i, j, k)*(qrn_and+qsn_and)
        tot_prec_anv = tot_prec_anv + (qrn_an+qsn_an)*mass(i, j, k)
        area_anv_prcd = area_anv_prcd + mass(i, j, k)*(cf_cond(i, j, k)*&
&         (qrn_an+qsn_an)+cf_con(i, j, k)*(qrn_and+qsn_and))
        area_anv_prc = area_anv_prc + cf_con(i, j, k)*(qrn_an+qsn_an)*&
&         mass(i, j, k)
        tot_prec_lsd = tot_prec_lsd + mass(i, j, k)*(qrn_lsd+qsn_lsd)
        tot_prec_ls = tot_prec_ls + (qrn_ls+qsn_ls)*mass(i, j, k)
        area_ls_prcd = area_ls_prcd + mass(i, j, k)*(cf_lsd(i, j, k)*(&
&         qrn_ls+qsn_ls)+cf_ls(i, j, k)*(qrn_lsd+qsn_lsd))
        area_ls_prc = area_ls_prc + cf_ls(i, j, k)*(qrn_ls+qsn_ls)*mass(&
&         i, j, k)
        IF (tot_prec_anv .GT. 0.0) THEN
          IF (area_anv_prc/tot_prec_anv .LT. 1.e-6) THEN
            area_anv_prc1 = 1.e-6
            area_anv_prc1d = 0.0_8
          ELSE
            area_anv_prc1d = (area_anv_prcd*tot_prec_anv-area_anv_prc*&
&             tot_prec_anvd)/tot_prec_anv**2
            area_anv_prc1 = area_anv_prc/tot_prec_anv
          END IF
        ELSE
          area_anv_prc1d = 0.0_8
        END IF
        IF (tot_prec_upd .GT. 0.0) THEN
          IF (area_upd_prc/tot_prec_upd .LT. 1.e-6) THEN
            area_upd_prc1 = 1.e-6
            area_upd_prc1d = 0.0_8
          ELSE
            area_upd_prc1d = (area_upd_prcd*tot_prec_upd-area_upd_prc*&
&             tot_prec_updd)/tot_prec_upd**2
            area_upd_prc1 = area_upd_prc/tot_prec_upd
          END IF
        ELSE
          area_upd_prc1d = 0.0_8
        END IF
        IF (tot_prec_ls .GT. 0.0) THEN
          IF (area_ls_prc/tot_prec_ls .LT. 1.e-6) THEN
            area_ls_prc1 = 1.e-6
            area_ls_prc1d = 0.0_8
          ELSE
            area_ls_prc1d = (area_ls_prcd*tot_prec_ls-area_ls_prc*&
&             tot_prec_lsd)/tot_prec_ls**2
            area_ls_prc1 = area_ls_prc/tot_prec_ls
          END IF
        ELSE
          area_ls_prc1d = 0.0_8
        END IF
        area_ls_prc1d = ls_beta*area_ls_prc1d
        area_ls_prc1 = ls_beta*area_ls_prc1
        area_upd_prc1d = cnv_beta*area_upd_prc1d
        area_upd_prc1 = cnv_beta*area_upd_prc1
        area_anv_prc1d = anv_beta*area_anv_prc1d
        area_anv_prc1 = anv_beta*area_anv_prc1
        IF (k .EQ. lm) THEN
! We have accumulated over the whole column
          IF (tot_prec_anv .GT. 0.0) THEN
            IF (area_anv_prc/tot_prec_anv .LT. 1.e-6) THEN
              area_anv_prc = 1.e-6
              area_anv_prcd = 0.0_8
            ELSE
              area_anv_prcd = (area_anv_prcd*tot_prec_anv-area_anv_prc*&
&               tot_prec_anvd)/tot_prec_anv**2
              area_anv_prc = area_anv_prc/tot_prec_anv
            END IF
          END IF
          IF (tot_prec_upd .GT. 0.0) THEN
            IF (area_upd_prc/tot_prec_upd .LT. 1.e-6) THEN
              area_upd_prc = 1.e-6
              area_upd_prcd = 0.0_8
            ELSE
              area_upd_prcd = (area_upd_prcd*tot_prec_upd-area_upd_prc*&
&               tot_prec_updd)/tot_prec_upd**2
              area_upd_prc = area_upd_prc/tot_prec_upd
            END IF
          END IF
          IF (tot_prec_ls .GT. 0.0) THEN
            IF (area_ls_prc/tot_prec_ls .LT. 1.e-6) THEN
              area_ls_prc = 1.e-6
              area_ls_prcd = 0.0_8
            ELSE
              area_ls_prcd = (area_ls_prcd*tot_prec_ls-area_ls_prc*&
&               tot_prec_lsd)/tot_prec_ls**2
              area_ls_prc = area_ls_prc/tot_prec_ls
            END IF
          END IF
          area_ls_prcd = ls_beta*area_ls_prcd
          area_ls_prc = ls_beta*area_ls_prc
          area_upd_prcd = cnv_beta*area_upd_prcd
          area_upd_prc = cnv_beta*area_upd_prc
          area_anv_prcd = anv_beta*area_anv_prcd
          area_anv_prc = anv_beta*area_anv_prc
        END IF
!Get micro-physical constants
        CALL cons_alhx_d(t(i, j, k), td(i, j, k), alhx3, alhx3d, &
&                  t_ice_max, t_ice_all, cons_alhs, cons_alhl)
        CALL cons_microphys_d(t(i, j, k), td(i, j, k), ph(i, j, k), qs(i&
&                       , j, k), qsd(i, j, k), aa, aad, bb, bbd, &
&                       cons_h2omw, cons_airmw, cons_rvap, alhx3, alhx3d&
&                      )
!Precip Scheme Expects Total Cloud Liquid
        qlt_tmpd = ql_lsd(i, j, k) + ql_cond(i, j, k)
        qlt_tmp = ql_ls(i, j, k) + ql_con(i, j, k)
        qit_tmpd = qi_lsd(i, j, k) + qi_cond(i, j, k)
        qit_tmp = qi_ls(i, j, k) + qi_con(i, j, k)
        prn_above_cu_oldd = prn_above_cu_newd
        prn_above_cu_old = prn_above_cu_new
        psn_above_cu_oldd = psn_above_cu_newd
        psn_above_cu_old = psn_above_cu_new
        evap_dd_cu_above_oldd = evap_dd_cu_above_newd
        evap_dd_cu_above_old = evap_dd_cu_above_new
        subl_dd_cu_above_oldd = subl_dd_cu_above_newd
        subl_dd_cu_above_old = subl_dd_cu_above_new
!Precip and Evap for Convection
        CALL precipandevap_d(k, ktop, lm, dt, frland(i, j), rhcrit, &
&                      qrn_cu_1d, qrn_cu_1dd, qsn_cu, qsn_cud, qlt_tmp, &
&                      qlt_tmpd, qit_tmp, t(i, j, k), td(i, j, k), q(i, &
&                      j, k), qd(i, j, k), mass(i, j, k), imass(i, j, k)&
&                      , ph(i, j, k), dzet(i, j, k), dzetd(i, j, k), &
&                      qddf3(i, j, k), qddf3d(i, j, k), aa, aad, bb, bbd&
&                      , area_upd_prc1, area_upd_prc1d, prn_above_cu_old&
&                      , prn_above_cu_oldd, prn_above_cu_new, &
&                      prn_above_cu_newd, psn_above_cu_old, &
&                      psn_above_cu_oldd, psn_above_cu_new, &
&                      psn_above_cu_newd, evap_dd_cu_above_old, &
&                      evap_dd_cu_above_oldd, evap_dd_cu_above_new, &
&                      evap_dd_cu_above_newd, subl_dd_cu_above_old, &
&                      subl_dd_cu_above_oldd, subl_dd_cu_above_new, &
&                      subl_dd_cu_above_newd, cnvenvfc, cnvddrfc, &
&                      cons_alhf, cons_alhs, cons_alhl, cons_cp, &
&                      cons_tice, cons_h2omw, cons_airmw, revap_off_p, &
&                      c_acc, c_ev_r, c_ev_s, rho_w, estblx)
        prn_above_an_oldd = prn_above_an_newd
        prn_above_an_old = prn_above_an_new
        psn_above_an_oldd = psn_above_an_newd
        psn_above_an_old = psn_above_an_new
        evap_dd_an_above_oldd = evap_dd_an_above_newd
        evap_dd_an_above_old = evap_dd_an_above_new
        subl_dd_an_above_oldd = subl_dd_an_above_newd
        subl_dd_an_above_old = subl_dd_an_above_new
!Precip and Evap for Anvil
        anvenvfc = 1.0
        CALL precipandevap_d(k, ktop, lm, dt, frland(i, j), rhcrit, &
&                      qrn_an, qrn_and, qsn_an, qsn_and, qlt_tmp, &
&                      qlt_tmpd, qit_tmp, t(i, j, k), td(i, j, k), q(i, &
&                      j, k), qd(i, j, k), mass(i, j, k), imass(i, j, k)&
&                      , ph(i, j, k), dzet(i, j, k), dzetd(i, j, k), &
&                      qddf3(i, j, k), qddf3d(i, j, k), aa, aad, bb, bbd&
&                      , area_anv_prc1, area_anv_prc1d, prn_above_an_old&
&                      , prn_above_an_oldd, prn_above_an_new, &
&                      prn_above_an_newd, psn_above_an_old, &
&                      psn_above_an_oldd, psn_above_an_new, &
&                      psn_above_an_newd, evap_dd_an_above_old, &
&                      evap_dd_an_above_oldd, evap_dd_an_above_new, &
&                      evap_dd_an_above_newd, subl_dd_an_above_old, &
&                      subl_dd_an_above_oldd, subl_dd_an_above_new, &
&                      subl_dd_an_above_newd, anvenvfc, anvddrfc, &
&                      cons_alhf, cons_alhs, cons_alhl, cons_cp, &
&                      cons_tice, cons_h2omw, cons_airmw, revap_off_p, &
&                      c_acc, c_ev_r, c_ev_s, rho_w, estblx)
        prn_above_ls_oldd = prn_above_ls_newd
        prn_above_ls_old = prn_above_ls_new
        psn_above_ls_oldd = psn_above_ls_newd
        psn_above_ls_old = psn_above_ls_new
        evap_dd_ls_above_oldd = evap_dd_ls_above_newd
        evap_dd_ls_above_old = evap_dd_ls_above_new
        subl_dd_ls_above_oldd = subl_dd_ls_above_newd
        subl_dd_ls_above_old = subl_dd_ls_above_new
!Precip and Evap for Large Scale
        lsenvfc = 1.0
        CALL precipandevap_d(k, ktop, lm, dt, frland(i, j), rhcrit, &
&                      qrn_ls, qrn_lsd, qsn_ls, qsn_lsd, qlt_tmp, &
&                      qlt_tmpd, qit_tmp, t(i, j, k), td(i, j, k), q(i, &
&                      j, k), qd(i, j, k), mass(i, j, k), imass(i, j, k)&
&                      , ph(i, j, k), dzet(i, j, k), dzetd(i, j, k), &
&                      qddf3(i, j, k), qddf3d(i, j, k), aa, aad, bb, bbd&
&                      , area_ls_prc1, area_ls_prc1d, prn_above_ls_old, &
&                      prn_above_ls_oldd, prn_above_ls_new, &
&                      prn_above_ls_newd, psn_above_ls_old, &
&                      psn_above_ls_oldd, psn_above_ls_new, &
&                      psn_above_ls_newd, evap_dd_ls_above_old, &
&                      evap_dd_ls_above_oldd, evap_dd_ls_above_new, &
&                      evap_dd_ls_above_newd, subl_dd_ls_above_old, &
&                      subl_dd_ls_above_oldd, subl_dd_ls_above_new, &
&                      subl_dd_ls_above_newd, lsenvfc, lsddrfc, &
&                      cons_alhf, cons_alhs, cons_alhl, cons_cp, &
&                      cons_tice, cons_h2omw, cons_airmw, revap_off_p, &
&                      c_acc, c_ev_r, c_ev_s, rho_w, estblx)
        IF (ql_ls(i, j, k) + ql_con(i, j, k) .GT. 0.00) THEN
          qt_tmpi_1d = -((ql_lsd(i, j, k)+ql_cond(i, j, k))/(ql_ls(i, j&
&           , k)+ql_con(i, j, k))**2)
          qt_tmpi_1 = 1./(ql_ls(i, j, k)+ql_con(i, j, k))
        ELSE
          qt_tmpi_1 = 0.0
          qt_tmpi_1d = 0.0_8
        END IF
        ql_lsd(i, j, k) = ql_lsd(i, j, k)*qlt_tmp*qt_tmpi_1 + ql_ls(i, j&
&         , k)*(qlt_tmpd*qt_tmpi_1+qlt_tmp*qt_tmpi_1d)
        ql_ls(i, j, k) = ql_ls(i, j, k)*qlt_tmp*qt_tmpi_1
        ql_cond(i, j, k) = ql_cond(i, j, k)*qlt_tmp*qt_tmpi_1 + ql_con(i&
&         , j, k)*(qlt_tmpd*qt_tmpi_1+qlt_tmp*qt_tmpi_1d)
        ql_con(i, j, k) = ql_con(i, j, k)*qlt_tmp*qt_tmpi_1
        IF (qi_ls(i, j, k) + qi_con(i, j, k) .GT. 0.00) THEN
          qt_tmpi_2d = -((qi_lsd(i, j, k)+qi_cond(i, j, k))/(qi_ls(i, j&
&           , k)+qi_con(i, j, k))**2)
          qt_tmpi_2 = 1./(qi_ls(i, j, k)+qi_con(i, j, k))
        ELSE
          qt_tmpi_2 = 0.0
          qt_tmpi_2d = 0.0_8
        END IF
        qi_lsd(i, j, k) = qi_lsd(i, j, k)*qit_tmp*qt_tmpi_2 + qi_ls(i, j&
&         , k)*(qit_tmpd*qt_tmpi_2+qit_tmp*qt_tmpi_2d)
        qi_ls(i, j, k) = qi_ls(i, j, k)*qit_tmp*qt_tmpi_2
        qi_cond(i, j, k) = qi_cond(i, j, k)*qit_tmp*qt_tmpi_2 + qi_con(i&
&         , j, k)*(qit_tmpd*qt_tmpi_2+qit_tmp*qt_tmpi_2d)
        qi_con(i, j, k) = qi_con(i, j, k)*qit_tmp*qt_tmpi_2

        !SINK FILTERING
        if (do_moist_physics == 1) then
           sinkfilt_ql  = 0.65
           sinkfilt_qi  = 0.65
           sinkfilt_cf  = 1.0
        elseif (do_moist_physics == 2) then
           sinkfilt_ql  = 0.9
           sinkfilt_qi  = 0.9
           sinkfilt_cf  = 1.0
        endif

        !Sink reduction on convective high altitude cloud lidquid ice cloud
        if (k < 50) then
            qi_lsd(i,j,k)  = sinkfilt_qi * qi_lsd(i,j,k)  +  (1.0-sinkfilt_qi) * qi_ls_p_presink
            qi_cond(i,j,k) = sinkfilt_qi * qi_cond(i,j,k) +  (1.0-sinkfilt_qi) * qi_con_p_presink
            qd(i,j,k)      = sinkfilt_qi * qd(i,j,k)      +  (1.0-sinkfilt_qi) * q_p_presink
        endif

        !Sink reduction on cloud liquid water
        if ( abs(k-62) .le. 2) then
           ql_lsd(i,j,k)  = sinkfilt_ql * ql_lsd(i,j,k)  +  (1.0-sinkfilt_ql) * ql_ls_p_presink
           ql_cond(i,j,k) = sinkfilt_ql * ql_cond(i,j,k) +  (1.0-sinkfilt_ql) * ql_con_p_presink
        endif

        !cf_cond(i,j,k) = SINKfilt_CF * cf_cond(i,j,k) +  (1.0-SINKfilt_CF) * cf_con_p_presink

        !TOTAL FILTERING
        totfilt_t  = 0.25
        td(i,j,k)      = totfilt_t  * td(i,j,k)      +  (1.0-totfilt_t) * t_p_preall

        if (do_moist_physics == 1) then
           totfilt_ql = 0.75
           totfilt_qi = 1.0
        elseif (do_moist_physics == 2) then
           totfilt_ql = 0.5
           totfilt_qi = 1.0!0.5
        endif

        ql_lsd(i,j,k)  = totfilt_ql * ql_lsd(i,j,k)  +  (1.0-totfilt_ql) * ql_ls_p_preall
        ql_cond(i,j,k) = totfilt_ql * ql_cond(i,j,k) +  (1.0-totfilt_ql) * ql_con_p_preall

        qi_lsd(i,j,k)  = totfilt_qi * qi_lsd(i,j,k)  +  (1.0-totfilt_qi) * qi_ls_p_preall
        qi_cond(i,j,k) = totfilt_qi * qi_cond(i,j,k) +  (1.0-totfilt_qi) * qi_con_p_preall

      END DO
    END DO
  END DO
!Clean up of excess relative humidity
  rhexcess = 1.1
  CALL dqsat_bac_d(dqsdt, dqsdtd, qs, qsd, t, td, ph, im, jm, lm, estblx&
&            , cons_h2omw, cons_airmw)
  dqsd = 0.0_8
  WHERE (q .GT. rhexcess*qs) 
    dqsd = ((qd-rhexcess*qsd)*(1.0+rhexcess*dqsdt*cons_alhl/cons_cp)-(q-&
&     rhexcess*qs)*rhexcess*cons_alhl*dqsdtd/cons_cp)/(1.0+rhexcess*&
&     dqsdt*cons_alhl/cons_cp)**2
    dqs = (q-rhexcess*qs)/(1.0+rhexcess*dqsdt*cons_alhl/cons_cp)
  ELSEWHERE
    dqsd = 0.0_8
    dqs = 0.0
  END WHERE
  qd = qd - dqsd
  q = q - dqs
  td = td + cons_alhl*dqsd/cons_cp
  t = t + cons_alhl/cons_cp*dqs
!Clean up Q<0
  DO j=1,jm
    DO i=1,im
!Total precipitable water
      tpwd = sum(dm(i, j, :)*qd(i, j, :))
      tpw = sum(q(i, j, :)*dm(i, j, :))
      negtpw = 0.
      negtpwd = 0.0_8
      DO l=1,lm
        IF (q(i, j, l) .LT. 0.0) THEN
          negtpwd = negtpwd + dm(i, j, l)*qd(i, j, l)
          negtpw = negtpw + q(i, j, l)*dm(i, j, l)
          qd(i, j, l) = 0.0_8
          q(i, j, l) = 0.0
        END IF
      END DO
      DO l=1,lm
        IF (q(i, j, l) .GE. 0.0) THEN
          qd(i, j, l) = qd(i, j, l)*(1.0+negtpw/(tpw-negtpw)) + q(i, j, &
&           l)*(negtpwd*(tpw-negtpw)-negtpw*(tpwd-negtpwd))/(tpw-negtpw)&
&           **2
          q(i, j, l) = q(i, j, l)*(1.0+negtpw/(tpw-negtpw))
        END IF
      END DO
    END DO
  END DO
!Convert temperature back to potential temperature
  thd = td/pih
  th = t/pih
END SUBROUTINE cloud_driver_d

!  Differentiation of cloud_tidy in forward (tangent) mode:
!   variations   of useful results: af qv qla qlc qia qic cf te
!   with respect to varying inputs: af qv qla qlc qia qic cf te
!SUBROUTINES
SUBROUTINE cloud_tidy_d(qv, qvd, te, ted, qlc, qlcd, qic, qicd, cf, cfd&
& , qla, qlad, qia, qiad, af, afd, cons_alhl, cons_alhs, cons_cp)
  IMPLICIT NONE
  REAL*8, INTENT(INOUT) :: te, qv, qlc, cf, qla, af, qic, qia
  REAL*8, INTENT(INOUT) :: ted, qvd, qlcd, cfd, qlad, afd, qicd, qiad
  REAL*8, INTENT(IN) :: cons_alhl, cons_alhs, cons_cp
!Fix if Anvil cloud fraction too small
  IF (af .LT. 1.e-5) THEN
    qvd = qvd + qlad + qiad
    qv = qv + qla + qia
    ted = ted - cons_alhl*qlad/cons_cp - cons_alhs*qiad/cons_cp
    te = te - cons_alhl/cons_cp*qla - cons_alhs/cons_cp*qia
    af = 0.
    qla = 0.
    qia = 0.
    afd = 0.0_8
    qlad = 0.0_8
    qiad = 0.0_8
  END IF
!Fix if LS cloud fraction too small
! if ( CF < 1.E-5 ) then
!         QV = QV + QLC + QIC
!         TE = TE - (CONS_ALHL/CONS_CP)*QLC - (CONS_ALHS/CONS_CP)*QIC
!         CF  = 0.
!         QLC = 0.
!         QIC = 0.
! end if
!LS LIQUID too small
  IF (qlc .LT. 1.e-8) THEN
    qvd = qvd + qlcd
    qv = qv + qlc
    ted = ted - cons_alhl*qlcd/cons_cp
    te = te - cons_alhl/cons_cp*qlc
    qlc = 0.
    qlcd = 0.0_8
  END IF
!LS ICE too small
  IF (qic .LT. 1.e-8) THEN
    qvd = qvd + qicd
    qv = qv + qic
    ted = ted - cons_alhs*qicd/cons_cp
    te = te - cons_alhs/cons_cp*qic
    qic = 0.
    qicd = 0.0_8
  END IF
!Anvil LIQUID too small
  IF (qla .LT. 1.e-8) THEN
    qvd = qvd + qlad
    qv = qv + qla
    ted = ted - cons_alhl*qlad/cons_cp
    te = te - cons_alhl/cons_cp*qla
    qla = 0.
    qlad = 0.0_8
  END IF
!Anvil ICE too small
  IF (qia .LT. 1.e-8) THEN
    qvd = qvd + qiad
    qv = qv + qia
    ted = ted - cons_alhs*qiad/cons_cp
    te = te - cons_alhs/cons_cp*qia
    qia = 0.
    qiad = 0.0_8
  END IF
!Fix ALL cloud quants if Anvil cloud LIQUID+ICE too small
  IF (qla + qia .LT. 1.e-8) THEN
    qvd = qvd + qlad + qiad
    qv = qv + qla + qia
    ted = ted - cons_alhl*qlad/cons_cp - cons_alhs*qiad/cons_cp
    te = te - cons_alhl/cons_cp*qla - cons_alhs/cons_cp*qia
    af = 0.
    qla = 0.
    qia = 0.
    afd = 0.0_8
    qlad = 0.0_8
    qiad = 0.0_8
  END IF
!Ditto if LS cloud LIQUID+ICE too small
  IF (qlc + qic .LT. 1.e-8) THEN
    qvd = qvd + qlcd + qicd
    qv = qv + qlc + qic
    ted = ted - cons_alhl*qlcd/cons_cp - cons_alhs*qicd/cons_cp
    te = te - cons_alhl/cons_cp*qlc - cons_alhs/cons_cp*qic
    cf = 0.
    qlc = 0.
    qic = 0.
    qlcd = 0.0_8
    qicd = 0.0_8
    cfd = 0.0_8
  END IF
END SUBROUTINE cloud_tidy_d

!  Differentiation of meltfreeze in forward (tangent) mode:
!   variations   of useful results: qi ql te
!   with respect to varying inputs: qi ql te
SUBROUTINE meltfreeze_d(dt, te, ted, ql, qld, qi, qid, t_ice_all, &
& t_ice_max, icefrpwr, cons_alhl, cons_alhs, cons_cp)
  IMPLICIT NONE
!Inputs
  REAL*8, INTENT(IN) :: dt, t_ice_all, t_ice_max
  INTEGER, INTENT(IN) :: icefrpwr
  REAL*8, INTENT(IN) :: cons_alhl, cons_alhs, cons_cp
!Prognostic
  REAL*8, INTENT(INOUT) :: te, ql, qi
  REAL*8, INTENT(INOUT) :: ted, qld, qid
!Locals
  REAL*8 :: fqi, dqil
  REAL*8 :: fqid, dqild
  REAL*8, PARAMETER :: taufrz=1000.
  INTRINSIC exp
  INTRINSIC max
  INTRINSIC min
  REAL*8 :: arg1
  REAL*8 :: arg1d
  fqi = 0.0
  dqil = 0.0
  CALL get_ice_fraction_d(te, ted, t_ice_all, t_ice_max, icefrpwr, fqi, &
&                   fqid)
!Freeze liquid
  IF (te .LE. t_ice_max) THEN
    arg1d = -(dt*fqid/taufrz)
    arg1 = -(dt*fqi/taufrz)
    dqild = qld*(1.0-exp(arg1)) - ql*arg1d*exp(arg1)
    dqil = ql*(1.0-exp(arg1))
  ELSE
    dqild = 0.0_8
  END IF
  IF (0. .LT. dqil) THEN
    dqil = dqil
  ELSE
    dqil = 0.
    dqild = 0.0_8
  END IF
  qid = qid + dqild
  qi = qi + dqil
  qld = qld - dqild
  ql = ql - dqil
  ted = ted + (cons_alhs-cons_alhl)*dqild/cons_cp
  te = te + (cons_alhs-cons_alhl)*dqil/cons_cp
  dqil = 0.
!Melt ice instantly above 0^C
  IF (te .GT. t_ice_max) THEN
    dqild = -qid
    dqil = -qi
  ELSE
    dqild = 0.0_8
  END IF
  IF (0. .GT. dqil) THEN
    dqil = dqil
  ELSE
    dqil = 0.
    dqild = 0.0_8
  END IF
  qid = qid + dqild
  qi = qi + dqil
  qld = qld - dqild
  ql = ql - dqil
  ted = ted + (cons_alhs-cons_alhl)*dqild/cons_cp
  te = te + (cons_alhs-cons_alhl)*dqil/cons_cp
END SUBROUTINE meltfreeze_d

!  Differentiation of convec_src in forward (tangent) mode:
!   variations   of useful results: af qv qla qia te
!   with respect to varying inputs: af qs qv qla qia dcf dmf te
SUBROUTINE convec_src_d(dt, mass, imass, te, ted, qv, qvd, dcf, dcfd, &
& dmf, dmfd, qla, qlad, qia, qiad, af, afd, qs, qsd, cons_alhs, &
& cons_alhl, cons_cp, t_ice_all, t_ice_max, icefrpwr)
  IMPLICIT NONE
!Inputs
  REAL*8, INTENT(IN) :: dt, t_ice_all, t_ice_max
  INTEGER, INTENT(IN) :: icefrpwr
  REAL*8, INTENT(IN) :: mass, imass, qs
  REAL*8, INTENT(IN) :: qsd
  REAL*8, INTENT(IN) :: dmf, dcf
  REAL*8, INTENT(IN) :: dmfd, dcfd
  REAL*8, INTENT(IN) :: cons_alhs, cons_alhl, cons_cp
!Prognostic
  REAL*8, INTENT(INOUT) :: te, qv
  REAL*8, INTENT(INOUT) :: ted, qvd
  REAL*8, INTENT(INOUT) :: qla, qia, af
  REAL*8, INTENT(INOUT) :: qlad, qiad, afd
!Locals
!Minimum allowed env RH
  REAL*8, PARAMETER :: minrhx=0.001
  REAL*8 :: tend, qvx, fqi
  REAL*8 :: tendd, fqid
  INTRINSIC min
!Namelist
!DT         - Timestep
!MASS       - Level Mass
!iMASS      - 1/Level Mass
!TE         - Temperature
!QV         - Specific Humidity
!DCF        - CNV_DQL from RAS
!DMF        - CNV_MFD from RAS
!QLA        - Convective cloud liquid water
!QIA        - Convective cloud liquid ice
!AF         - Convective cloud fraction
!QS         - Qsat
!Zero out locals
  tend = 0.0
  qvx = 0.0
  fqi = 0.0
!Addition of condensate from RAS
  tendd = imass*dcfd
  tend = dcf*imass
  CALL get_ice_fraction_d(te, ted, t_ice_all, t_ice_max, icefrpwr, fqi, &
&                   fqid)
  qlad = qlad + dt*((1.0-fqi)*tendd-fqid*tend)
  qla = qla + (1.0-fqi)*tend*dt
  qiad = qiad + dt*(fqid*tend+fqi*tendd)
  qia = qia + fqi*tend*dt
!Convective condensation has never frozen so latent heat of fusion
  ted = ted + (cons_alhs-cons_alhl)*dt*(fqid*tend+fqi*tendd)/cons_cp
  te = te + (cons_alhs-cons_alhl)*fqi*tend*dt/cons_cp
!Compute Tiedtke-style anvil fraction
  tendd = imass*dmfd
  tend = dmf*imass
  afd = afd + dt*tendd
  af = af + tend*dt
  IF (af .GT. 0.99) THEN
    af = 0.99
    afd = 0.0_8
  ELSE
    af = af
  END IF
! Check for funny (tiny, negative) external QV, resulting from assumed QV=QSAT within anvil.
! Simply constrain AF assume condensate just gets "packed" in     
  IF (af .LT. 1.0) THEN
    qvx = (qv-qs*af)/(1.-af)
  ELSE
    qvx = qs
  END IF
!If saturated over critial value and there is Anvil
  IF (qvx - minrhx*qs .LT. 0.0 .AND. af .GT. 0.) THEN
    afd = ((qvd-minrhx*qsd)*qs*(1.0-minrhx)-(qv-minrhx*qs)*(1.0-minrhx)*&
&     qsd)/(qs*(1.0-minrhx))**2
    af = (qv-minrhx*qs)/(qs*(1.0-minrhx))
  END IF
  IF (af .LT. 0.) THEN
! If still cant make suitable env RH then destroy anvil
    af = 0.0
    qvd = qvd + qlad + qiad
    qv = qv + qla + qia
    ted = ted - (cons_alhl*qlad+cons_alhs*qiad)/cons_cp
    te = te - (cons_alhl*qla+cons_alhs*qia)/cons_cp
    qla = 0.0
    qia = 0.0
    afd = 0.0_8
    qlad = 0.0_8
    qiad = 0.0_8
  END IF
END SUBROUTINE convec_src_d

!  Differentiation of ls_cloud in forward (tangent) mode:
!   variations   of useful results: qai qal af qv qci qcl cf te
!   with respect to varying inputs: qai qal af qv qci qcl cf te
SUBROUTINE ls_cloud_d(dt, alpha, pdfshape, pl, te, ted, qv, qvd, qcl, &
& qcld, qal, qald, qci, qcid, qai, qaid, cf, cfd, af, afd, cons_alhl, &
& cons_alhf, cons_alhs, cons_cp, cons_h2omw, cons_airmw, t_ice_all, &
& t_ice_max, icefrpwr, estblx, cloud_pertmod, dmp)
  IMPLICIT NONE
!Inputs
  REAL*8, INTENT(IN) :: dt, alpha, pl, t_ice_all, t_ice_max
  INTEGER, INTENT(IN) :: pdfshape, cloud_pertmod, dmp
  INTEGER, INTENT(IN) :: icefrpwr
  REAL*8, INTENT(IN) :: cons_alhl, cons_alhf, cons_alhs, cons_cp, &
& cons_h2omw, cons_airmw
  REAL*8, INTENT(IN) :: estblx(:)
!Prognostic
  REAL*8, INTENT(INOUT) :: te, qv, qcl, qci, qal, qai, cf, af
  REAL*8, INTENT(INOUT) :: ted, qvd, qcld, qcid, qald, qaid, cfd, afd
!Locals
  INTEGER :: n
  REAL*8 :: qco, cfo, qao, qt, qmx, qmn, dq
  REAL*8 :: qcod, cfod, qaod, qtd
  REAL*8 :: teo, qsx, dqsx, qs, dqs, tmparr
  REAL*8 :: teod, qsxd, dqsxd, dqsd, tmparrd
  REAL*8 :: qcx, qvx, cfx, qax, qc, qa, fqi, fqi_a, dqai, dqal, dqci, &
& dqcl
  REAL*8 :: qcxd, qvxd, cfxd, qaxd, qcd, qad, fqid, dqaid, dqald, dqcid&
& , dqcld
  REAL*8 :: ten, qsp, cfp, qvp, qcp
  REAL*8 :: tend, qcpd
  REAL*8 :: tep, qsn, cfn, qvn, qcn
  REAL*8 :: tepd, qsnd, cfnd, qcnd
  REAL*8 :: alhx, sigmaqt1, sigmaqt2
  REAL*8 :: alhxd, sigmaqt1d, sigmaqt2d
  REAL*8, DIMENSION(1) :: dqsx1, qsx1, teo1, pl1
  INTRINSIC max
!Namelist
!DT      - Timestep 
!ALPHA   - PDF half width
!PL      - Pressure (hPa)
!TE      - Temperature
!QV      - Specific humidity
!QCl     - Convective cloud liquid water
!QAl     - Large scale cloud liquid water
!QCi     - Convective cloud liquid ice
!QAi     - Large scale cloud liquid ice
!CF      - Large scale cloud fraction
!AF      - Convective cloud fraction
  qcd = qcld + qcid
  qc = qcl + qci
  qad = qald + qaid
  qa = qal + qai
  IF (qa .GT. 0.0) THEN
    fqi_a = qai/qa
  ELSE
    fqi_a = 0.0
  END IF
  teod = ted
  teo = te
  CALL dqsats_bac_d(dqsx, dqsxd, qsx, qsxd, teo, teod, pl, estblx, &
&             cons_h2omw, cons_airmw)
  IF (af .LT. 1.0) THEN
    IF (dmp .EQ. 1) THEN
      IF (1. - af .GT. 0.02) THEN
        tmparrd = -((-afd)/(1.-af)**2)
        tmparr = 1./(1.-af)
      ELSE
        tmparrd = -((-afd)/(0.02)**2)
        tmparr = 1./(1.-af)
      END IF
    ELSE IF (dmp .EQ. 2) THEN
      tmparrd = -((-afd)/(1.-af)**2)
      tmparr = 1./(1.-af)
    ELSE
      tmparrd = 0.0_8
    END IF
  ELSE
    tmparr = 0.0
    tmparrd = 0.0_8
  END IF
  cfxd = cfd*tmparr + cf*tmparrd
  cfx = cf*tmparr
  qcxd = qcd*tmparr + qc*tmparrd
  qcx = qc*tmparr
  qvxd = (qvd-qsxd*af-qsx*afd)*tmparr + (qv-qsx*af)*tmparrd
  qvx = (qv-qsx*af)*tmparr
  IF (af .GE. 1.0) THEN
    qvxd = 1.e-4*qsxd
    qvx = qsx*1.e-4
  END IF
  IF (af .GT. 0.) THEN
    qaxd = (qad*af-qa*afd)/af**2
    qax = qa/af
  ELSE
    qax = 0.
    qaxd = 0.0_8
  END IF
  qtd = qcxd + qvxd
  qt = qcx + qvx
  tep = teo
  qsn = qsx
  tend = teod
  ten = teo
  cfnd = cfxd
  cfn = cfx
  qvn = qvx
  qcnd = qcxd
  qcn = qcx
  dqs = dqsx
!Begin iteration
!do n=1,4
  n = 1
  qsp = qsn
  qvp = qvn
  qcpd = qcnd
  qcp = qcn
  cfp = cfn
!Dont call again as not looping
  dqsd = dqsxd
  dqs = dqsx
  qsnd = qsxd
  qsn = qsx
!call DQSATs_BAC(DQS, QSn, TEn, PL, ESTBLX, CONS_H2OMW, CONS_AIRMW)
  tepd = tend
  tep = ten
  CALL get_ice_fraction_d(tep, tepd, t_ice_all, t_ice_max, icefrpwr, fqi&
&                   , fqid)
  sigmaqt1d = alpha*qsnd
  sigmaqt1 = alpha*qsn
  sigmaqt2d = alpha*qsnd
  sigmaqt2 = alpha*qsn
  IF (pdfshape .EQ. 2) THEN
!For triangular, symmetric: sigmaqt1 = sigmaqt2 = alpha*qsn (alpha is half width)
!For triangular, skewed r : sigmaqt1 < sigmaqt2
    sigmaqt1d = alpha*qsnd
    sigmaqt1 = alpha*qsn
    sigmaqt2d = alpha*qsnd
    sigmaqt2 = alpha*qsn
  END IF
!Compute cloud fraction
  IF (cloud_pertmod .EQ. 0) THEN
    CALL pdffrac_d(1, qt, qtd, sigmaqt1, sigmaqt1d, sigmaqt2, sigmaqt2d&
&            , qsn, qsnd, cfn, cfnd)
  ELSE IF (cloud_pertmod .EQ. 1) THEN
    CALL pdffrac_d(4, qt, qtd, sigmaqt1, sigmaqt1d, sigmaqt2, sigmaqt2d&
&            , qsn, qsnd, cfn, cfnd)
  END IF
!Compute cloud condensate
  CALL pdfcondensate_d(pdfshape, qt, qtd, sigmaqt1, sigmaqt1d, sigmaqt2&
&                , sigmaqt2d, qsn, qsnd, qcn, qcnd)
!Adjustments to anvil condensate due to the assumption of a stationary TOTAL 
!water PDF subject to a varying QSAT value during the iteration
  IF (af .GT. 0.) THEN
! + QSx - QS 
    qaod = qaxd
    qao = qax
  ELSE
    qao = 0.
    qaod = 0.0_8
  END IF
  alhxd = cons_alhs*fqid - cons_alhl*fqid
  alhx = (1.0-fqi)*cons_alhl + fqi*cons_alhs
  IF (pdfshape .EQ. 1) THEN
    qcnd = qcpd + ((qcnd-qcpd)*(1.-(cfn*(alpha-1.)-qcn/qsn)*dqs*alhx/&
&     cons_cp)+(qcn-qcp)*(((alpha-1.)*cfnd-(qcnd*qsn-qcn*qsnd)/qsn**2)*&
&     dqs*alhx+(cfn*(alpha-1.)-qcn/qsn)*(dqsd*alhx+dqs*alhxd))/cons_cp)/&
&     (1.-(cfn*(alpha-1.)-qcn/qsn)*dqs*alhx/cons_cp)**2
    qcn = qcp + (qcn-qcp)/(1.-(cfn*(alpha-1.)-qcn/qsn)*dqs*alhx/cons_cp)
  ELSE IF (pdfshape .EQ. 2) THEN
!This next line needs correcting - need proper d(del qc)/dT derivative for triangular
!for now, just use relaxation of 1/2.
    IF (n .NE. 4) THEN
      qcnd = qcpd + 0.5*(qcnd-qcpd)
      qcn = qcp + (qcn-qcp)*0.5
    END IF
  END IF
  qvn = qvp - (qcn-qcp)
  ten = tep + (1.0-fqi)*(cons_alhl/cons_cp)*((qcn-qcp)*(1.-af)+(qao-qax)&
&   *af) + fqi*(cons_alhs/cons_cp)*((qcn-qcp)*(1.-af)+(qao-qax)*af)
!enddo ! qsat iteration
  cfod = cfnd
  cfo = cfn
  cf = cfn
  qcod = qcnd
  qco = qcn
  teo = ten
! Update prognostic variables. QCo, QAo become updated grid means.
  IF (af .LT. 1.0) THEN
    cfd = cfod*(1.-af) - cfo*afd
    cf = cfo*(1.-af)
    qcod = qcod*(1.-af) - qco*afd
    qco = qco*(1.-af)
    qaod = qaod*af + qao*afd
    qao = qao*af
  ELSE
!Grid box filled with anvil
! Special case AF=1, i.e., box filled with anvil. Note: no guarantee QV_box > QS_box
! Remove any other cloud
    cf = 0.
! Add any LS condensate to anvil type
    qaod = qad + qcd
    qao = qa + qc
! Remove same from LS   
    qco = 0.
! Total water
    qtd = qaod + qvd
    qt = qao + qv
    IF (qt - qsx .LT. 0.) THEN
      qao = 0.
      cfd = 0.0_8
      qaod = 0.0_8
      qcod = 0.0_8
    ELSE
      qaod = qtd - qsxd
      qao = qt - qsx
      cfd = 0.0_8
      qcod = 0.0_8
    END IF
  END IF
!Partition new condensate into ice and liquid taking care to keep both >=0 separately. 
!New condensate can be less than old, so Delta can be < 0.
  qcxd = qcod - qcd
  qcx = qco - qc
  dqcld = (1.0-fqi)*qcxd - fqid*qcx
  dqcl = (1.0-fqi)*qcx
  dqcid = fqid*qcx + fqi*qcxd
  dqci = fqi*qcx
!Large Scale Partition
  IF (qcl + dqcl .LT. 0.) THEN
    dqcid = dqcid + qcld + dqcld
    dqci = dqci + (qcl+dqcl)
!== dQCl - (QCl+dQCl)
    dqcld = -qcld
    dqcl = -qcl
  END IF
  IF (qci + dqci .LT. 0.) THEN
    dqcld = dqcld + qcid + dqcid
    dqcl = dqcl + (qci+dqci)
!== dQCi - (QCi+dQCi)
    dqcid = -qcid
    dqci = -qci
  END IF
  qaxd = qaod - qad
  qax = qao - qa
! (1.0-fQi)*QAx
  dqald = qaxd
  dqal = qax
!  fQi  * QAx
  dqai = 0.
!Convective partition
  IF (qal + dqal .LT. 0.) THEN
    dqaid = qald + dqald
    dqai = dqai + (qal+dqal)
    dqald = -qald
    dqal = -qal
  ELSE
    dqaid = 0.0_8
  END IF
  IF (qai + dqai .LT. 0.) THEN
    dqald = dqald + qaid + dqaid
    dqal = dqal + (qai+dqai)
    dqaid = -qaid
    dqai = -qai
  END IF
! Clean-up cloud if fractions are too small
  IF (af .LT. 1.e-5) THEN
    dqaid = -qaid
    dqai = -qai
    dqald = -qald
    dqal = -qal
  END IF
  IF (cf .LT. 1.e-5) THEN
    dqcid = -qcid
    dqci = -qci
    dqcld = -qcld
    dqcl = -qcl
  END IF
  qaid = qaid + dqaid
  qai = qai + dqai
  qald = qald + dqald
  qal = qal + dqal
  qcid = qcid + dqcid
  qci = qci + dqci
  qcld = qcld + dqcld
  qcl = qcl + dqcl
!Update specific humidity
  qvd = qvd - dqaid - dqcid - dqald - dqcld
  qv = qv - (dqai+dqci+dqal+dqcl)
!Update temperature
  ted = ted + (cons_alhl*(dqaid+dqcid+dqald+dqcld)+cons_alhf*(dqaid+&
&   dqcid))/cons_cp
  te = te + (cons_alhl*(dqai+dqci+dqal+dqcl)+cons_alhf*(dqai+dqci))/&
&   cons_cp
!Take care of situations where QS moves past QA during QSAT iteration (QAo negative). 
!"Evaporate" offending QA
  IF (qao .LE. 0.) THEN
    qvd = qvd + qaid + qald
    qv = qv + qai + qal
    ted = ted - cons_alhs*qaid/cons_cp - cons_alhl*qald/cons_cp
    te = te - cons_alhs/cons_cp*qai - cons_alhl/cons_cp*qal
    qai = 0.
    qal = 0.
    af = 0.
    qaid = 0.0_8
    qald = 0.0_8
    afd = 0.0_8
  END IF
END SUBROUTINE ls_cloud_d

!  Differentiation of pdffrac in forward (tangent) mode:
!   variations   of useful results: clfrac
!   with respect to varying inputs: qtmean sigmaqt1 sigmaqt2 qstar
!                clfrac
SUBROUTINE pdffrac_d(flag, qtmean, qtmeand, sigmaqt1, sigmaqt1d, &
& sigmaqt2, sigmaqt2d, qstar, qstard, clfrac, clfracd)
  IMPLICIT NONE
!Regularization
!clfracd = 0.2*clfracd
!Inputs
  INTEGER, INTENT(IN) :: flag
  REAL*8, INTENT(IN) :: qtmean, sigmaqt1, sigmaqt2, qstar
  REAL*8, INTENT(IN) :: qtmeand, sigmaqt1d, sigmaqt2d, qstard
!Prognostic
  REAL*8, INTENT(INOUT) :: clfrac
  REAL*8, INTENT(INOUT) :: clfracd
!LOCALS
  REAL*8 :: qtmode, qtmin, qtmax
  REAL*8 :: qtmoded, qtmind, qtmaxd
  REAL*8 :: rh, rhd, q1, q2
  REAL*8 :: rhd0
  INTRINSIC min
  INTRINSIC tanh
  REAL*8 :: min1
  REAL*8 :: min1d
  IF (flag .EQ. 1) THEN
!Tophat PDF
    IF (qtmean + sigmaqt1 .LT. qstar) THEN
      clfrac = 0.
      clfracd = 0.0_8
    ELSE IF (sigmaqt1 .GT. 0.) THEN
      IF (qtmean + sigmaqt1 - qstar .GT. 2.*sigmaqt1) THEN
        min1d = 2.*sigmaqt1d
        min1 = 2.*sigmaqt1
      ELSE
        min1d = qtmeand + sigmaqt1d - qstard
        min1 = qtmean + sigmaqt1 - qstar
      END IF
      clfracd = (min1d*2.*sigmaqt1-min1*2.*sigmaqt1d)/(2.*sigmaqt1)**2
      clfrac = min1/(2.*sigmaqt1)
    ELSE
      clfrac = 1.
      clfracd = 0.0_8
    END IF
  ELSE IF (flag .EQ. 2) THEN
!Triangular PDF
    qtmoded = qtmeand + (sigmaqt1d-sigmaqt2d)/3.
    qtmode = qtmean + (sigmaqt1-sigmaqt2)/3.
    IF (qtmode - sigmaqt1 .GT. 0.) THEN
      qtmin = 0.
      qtmind = 0.0_8
    ELSE
      qtmind = qtmoded - sigmaqt1d
      qtmin = qtmode - sigmaqt1
    END IF
    qtmaxd = qtmoded + sigmaqt2d
    qtmax = qtmode + sigmaqt2
    IF (qtmax .LT. qstar) THEN
      clfrac = 0.
      clfracd = 0.0_8
    ELSE IF (qtmode .LE. qstar .AND. qstar .LT. qtmax) THEN
      clfracd = (((qtmaxd-qstard)*(qtmax-qstar)+(qtmax-qstar)*(qtmaxd-&
&       qstard))*(qtmax-qtmin)*(qtmax-qtmode)-(qtmax-qstar)**2*((qtmaxd-&
&       qtmind)*(qtmax-qtmode)+(qtmax-qtmin)*(qtmaxd-qtmoded)))/((qtmax-&
&       qtmin)*(qtmax-qtmode))**2
      clfrac = (qtmax-qstar)*(qtmax-qstar)/((qtmax-qtmin)*(qtmax-qtmode)&
&       )
    ELSE IF (qtmin .LE. qstar .AND. qstar .LT. qtmode) THEN
      clfracd = -((((qstard-qtmind)*(qstar-qtmin)+(qstar-qtmin)*(qstard-&
&       qtmind))*(qtmax-qtmin)*(qtmode-qtmin)-(qstar-qtmin)**2*((qtmaxd-&
&       qtmind)*(qtmode-qtmin)+(qtmax-qtmin)*(qtmoded-qtmind)))/((qtmax-&
&       qtmin)*(qtmode-qtmin))**2)
      clfrac = 1. - (qstar-qtmin)*(qstar-qtmin)/((qtmax-qtmin)*(qtmode-&
&       qtmin))
    ELSE IF (qstar .LE. qtmin) THEN
      clfrac = 1.
      clfracd = 0.0_8
    END IF
  ELSE IF (flag .EQ. 3) THEN

    !Tophat PDF for the reference part
    IF (qtmean + sigmaqt1 .LT. qstar) THEN
      clfrac = 0.
    ELSE IF (sigmaqt1 .GT. 0.) THEN
      IF (qtmean + sigmaqt1 - qstar .GT. 2.*sigmaqt1) THEN
        min1d = 2.*sigmaqt1d
        min1 = 2.*sigmaqt1
      ELSE
        min1d = qtmeand + sigmaqt1d - qstard
        min1 = qtmean + sigmaqt1 - qstar
      END IF
      clfrac = min1/(2.*sigmaqt1)
    ELSE
      clfrac = 1.
    END IF

    !TANH function for the perturabtions
    rhd0 = (qtmeand*qstar-qtmean*qstard)/qstar**2
    rh = qtmean/qstar
    q1 = 22.556
    clfracd = 0.5*q1*rhd0*(1.0-tanh(q1*(rh-1.0))**2)

    ! (REGULARIZATION) * (GRADIENT) * (PERTURBATION)
    clfracd = (0.66*( cosh(10*(rh-1.0))**-2)) * clfracd

  ELSE IF (flag .EQ. 4) THEN

    !Tophat PDF for the reference part
    IF (qtmean + sigmaqt1 .LT. qstar) THEN
      clfrac = 0.
    ELSE IF (sigmaqt1 .GT. 0.) THEN
      IF (qtmean + sigmaqt1 - qstar .GT. 2.*sigmaqt1) THEN
        min1d = 2.*sigmaqt1d
        min1 = 2.*sigmaqt1
      ELSE
        min1d = qtmeand + sigmaqt1d - qstard
        min1 = qtmean + sigmaqt1 - qstar
      END IF
      clfrac = min1/(2.*sigmaqt1)
    ELSE
      clfrac = 1.
    END IF

    !Linear for the perturbation part
    rhd0 = (qtmeand*qstar-qtmean*qstard)/qstar**2
    rh = qtmean/qstar
    q1 = 0.9335
    q2 = 1.0665
    IF (rh .LT. q1) THEN
      clfracd = 0.0_8
    ELSE IF (rh .GE. q1 .AND. rh .LT. q2) THEN
      clfracd = rhd0/((q2/q1-1)*q1)
    ELSE
      clfracd = 0.0_8
    END IF

    !Regularization
    clfracd = clfracd * 0.2

  END IF

END SUBROUTINE pdffrac_d

!  Differentiation of pdfcondensate in forward (tangent) mode:
!   variations   of useful results: condensate4
!   with respect to varying inputs: qtmean4 qstar4 sigmaqt14 sigmaqt24
SUBROUTINE pdfcondensate_d(flag, qtmean4, qtmean4d, sigmaqt14, &
& sigmaqt14d, sigmaqt24, sigmaqt24d, qstar4, qstar4d, condensate4, &
& condensate4d)
  IMPLICIT NONE
!Inputs
  INTEGER, INTENT(IN) :: flag
  REAL*8, INTENT(IN) :: qtmean4, sigmaqt14, sigmaqt24, qstar4
  REAL*8, INTENT(IN) :: qtmean4d, sigmaqt14d, sigmaqt24d, qstar4d
!Prognostic
  REAL*8, INTENT(INOUT) :: condensate4
  REAL*8, INTENT(INOUT) :: condensate4d
!Locals
  REAL*8 :: qtmode, qtmin, qtmax, consta, constb, cloudf
  REAL*8 :: qtmoded, qtmind, qtmaxd, constad, constbd, cloudfd
  REAL*8 :: term1, term2, term3
  REAL*8 :: term1d, term2d, term3d
  REAL*8 :: qtmean, sigmaqt1, sigmaqt2, qstar, condensate
  REAL*8 :: qtmeand, sigmaqt1d, sigmaqt2d, qstard, condensated
  INTRINSIC min
  REAL*8 :: min1
  REAL*8 :: min1d
  qtmeand = qtmean4d
  qtmean = qtmean4
  sigmaqt1d = sigmaqt14d
  sigmaqt1 = sigmaqt14
  sigmaqt2d = sigmaqt24d
  sigmaqt2 = sigmaqt24
  qstard = qstar4d
  qstar = qstar4
  IF (flag .EQ. 1) THEN
    IF (qtmean + sigmaqt1 .LT. qstar) THEN
      condensate = 0.d0
      condensated = 0.0_8
    ELSE IF (qstar .GT. qtmean - sigmaqt1) THEN
      IF (sigmaqt1 .GT. 0.d0) THEN
        IF (qtmean + sigmaqt1 - qstar .GT. 2.d0*sigmaqt1) THEN
          min1d = 2.d0*sigmaqt1d
          min1 = 2.d0*sigmaqt1
        ELSE
          min1d = qtmeand + sigmaqt1d - qstard
          min1 = qtmean + sigmaqt1 - qstar
        END IF
        condensated = (2*min1*min1d*4.d0*sigmaqt1-min1**2*4.d0*sigmaqt1d&
&         )/(4.d0*sigmaqt1)**2
        condensate = min1**2/(4.d0*sigmaqt1)
      ELSE
        condensated = qtmeand - qstard
        condensate = qtmean - qstar
      END IF
    ELSE
      condensated = qtmeand - qstard
      condensate = qtmean - qstar
    END IF
  ELSE IF (flag .EQ. 2) THEN
    qtmoded = qtmeand + (sigmaqt1d-sigmaqt2d)/3.d0
    qtmode = qtmean + (sigmaqt1-sigmaqt2)/3.d0
    IF (qtmode - sigmaqt1 .GT. 0.d0) THEN
      qtmin = 0.d0
      qtmind = 0.0_8
    ELSE
      qtmind = qtmoded - sigmaqt1d
      qtmin = qtmode - sigmaqt1
    END IF
    qtmaxd = qtmoded + sigmaqt2d
    qtmax = qtmode + sigmaqt2
    IF (qtmax .LT. qstar) THEN
      condensate = 0.d0
      condensated = 0.0_8
    ELSE IF (qtmode .LE. qstar .AND. qstar .LT. qtmax) THEN
      constbd = -(2.d0*((qtmaxd-qtmind)*(qtmax-qtmode)+(qtmax-qtmin)*(&
&       qtmaxd-qtmoded))/((qtmax-qtmin)*(qtmax-qtmode))**2)
      constb = 2.d0/((qtmax-qtmin)*(qtmax-qtmode))
      cloudfd = (((qtmaxd-qstard)*(qtmax-qstar)+(qtmax-qstar)*(qtmaxd-&
&       qstard))*(qtmax-qtmin)*(qtmax-qtmode)-(qtmax-qstar)**2*((qtmaxd-&
&       qtmind)*(qtmax-qtmode)+(qtmax-qtmin)*(qtmaxd-qtmoded)))/((qtmax-&
&       qtmin)*(qtmax-qtmode))**2
      cloudf = (qtmax-qstar)*(qtmax-qstar)/((qtmax-qtmin)*(qtmax-qtmode)&
&       )
      term1d = ((qstard*qstar+qstar*qstard)*qstar+qstar**2*qstard)/3.d0
      term1 = qstar*qstar*qstar/3.d0
      term2d = ((qtmaxd*qstar+qtmax*qstard)*qstar+qtmax*qstar*qstard)/&
&       2.d0
      term2 = qtmax*qstar*qstar/2.d0
      term3d = ((qtmaxd*qtmax+qtmax*qtmaxd)*qtmax+qtmax**2*qtmaxd)/6.d0
      term3 = qtmax*qtmax*qtmax/6.d0
      condensated = constbd*(term1-term2+term3) + constb*(term1d-term2d+&
&       term3d) - qstard*cloudf - qstar*cloudfd
      condensate = constb*(term1-term2+term3) - qstar*cloudf
    ELSE IF (qtmin .LE. qstar .AND. qstar .LT. qtmode) THEN
      constad = -(2.d0*((qtmaxd-qtmind)*(qtmode-qtmin)+(qtmax-qtmin)*(&
&       qtmoded-qtmind))/((qtmax-qtmin)*(qtmode-qtmin))**2)
      consta = 2.d0/((qtmax-qtmin)*(qtmode-qtmin))
      cloudfd = -((((qstard-qtmind)*(qstar-qtmin)+(qstar-qtmin)*(qstard-&
&       qtmind))*(qtmax-qtmin)*(qtmode-qtmin)-(qstar-qtmin)**2*((qtmaxd-&
&       qtmind)*(qtmode-qtmin)+(qtmax-qtmin)*(qtmoded-qtmind)))/((qtmax-&
&       qtmin)*(qtmode-qtmin))**2)
      cloudf = 1.d0 - (qstar-qtmin)*(qstar-qtmin)/((qtmax-qtmin)*(qtmode&
&       -qtmin))
      term1d = ((qstard*qstar+qstar*qstard)*qstar+qstar**2*qstard)/3.d0
      term1 = qstar*qstar*qstar/3.d0
      term2d = ((qtmind*qstar+qtmin*qstard)*qstar+qtmin*qstar*qstard)/&
&       2.d0
      term2 = qtmin*qstar*qstar/2.d0
      term3d = ((qtmind*qtmin+qtmin*qtmind)*qtmin+qtmin**2*qtmind)/6.d0
      term3 = qtmin*qtmin*qtmin/6.d0
      condensated = qtmeand - constad*(term1-term2+term3) - consta*(&
&       term1d-term2d+term3d) - qstard*cloudf - qstar*cloudfd
      condensate = qtmean - consta*(term1-term2+term3) - qstar*cloudf
    ELSE IF (qstar .LE. qtmin) THEN
      condensated = qtmeand - qstard
      condensate = qtmean - qstar
    ELSE
      condensated = 0.0_8
    END IF
  ELSE IF (flag .EQ. 3) THEN
!Reference part done normally
!IF (qtmean + sigmaqt1 .LT. qstar) THEN
!  condensate = 0.d0
!ELSE IF (qstar .GT. qtmean - sigmaqt1) THEN
!  IF (sigmaqt1 .GT. 0.d0) THEN
!    IF (qtmean + sigmaqt1 - qstar .GT. 2.d0*sigmaqt1) THEN
!      min1 = 2.d0*sigmaqt1
!    ELSE
!      min1 = qtmean + sigmaqt1 - qstar
!    END IF
!    condensate = min1**2/(4.d0*sigmaqt1)
!  ELSE
!    condensate = qtmean - qstar
!  END IF
!ELSE
!  condensate = qtmean - qstar
!END IF
!Perturbation part from linear
    IF (qtmean - qstar .GT. -0.5e-3) THEN
      condensated = qtmeand - qstard
      condensate = qtmean - qstar
    ELSE
      condensate = 0.0
      condensated = 0.0_8
    END IF
  ELSE
    condensated = 0.0_8
  END IF
  condensate4d = condensated
  condensate4 = condensate
END SUBROUTINE pdfcondensate_d

!  Differentiation of evap_cnv in forward (tangent) mode:
!   variations   of useful results: f ql qv te
!   with respect to varying inputs: f qi ql qs qv te
SUBROUTINE evap_cnv_d(dt, rhcr, pl, te, ted, qv, qvd, ql, qld, qi, qid, &
& f, fd, xf, qs, qsd, rho_w, cld_evp_eff, cons_h2omw, cons_airmw, &
& cons_alhl, cons_rvap, cons_rgas, cons_pi, cons_cp)
  IMPLICIT NONE
!Inputs
  REAL*8, INTENT(IN) :: dt, rhcr, pl, xf, qs, rho_w, cld_evp_eff
  REAL*8, INTENT(IN) :: qsd
  REAL*8, INTENT(IN) :: cons_h2omw, cons_airmw, cons_alhl, cons_rvap, &
& cons_rgas, cons_pi, cons_cp
!Prognostics
  REAL*8, INTENT(INOUT) :: te, qv, ql, qi, f
  REAL*8, INTENT(INOUT) :: ted, qvd, qld, qid, fd
!Locals
  REAL*8 :: es, radius, k1, k2, teff, qcm, evap, rhx, qc, a_eff, epsilon
  REAL*8 :: esd, radiusd, k1d, k2d, teffd, qcmd, evapd, rhxd, qcd
  REAL*8, PARAMETER :: k_cond=2.4e-2
  REAL*8, PARAMETER :: diffu=2.2e-5
  REAL*8, PARAMETER :: nn=50.*1.0e6
  INTRINSIC min
  epsilon = cons_h2omw/cons_airmw
  a_eff = cld_evp_eff
!EVAPORATION OF CLOUD WATER.
! (100 <-^ convert from mbar to Pa)
  esd = (100.*pl*qsd*(epsilon+(1.0-epsilon)*qs)-100.*pl*qs*(1.0-epsilon)&
&   *qsd)/(epsilon+(1.0-epsilon)*qs)**2
  es = 100.*pl*qs/(epsilon+(1.0-epsilon)*qs)
  IF (qv/qs .GT. 1.00) THEN
    rhx = 1.00
    rhxd = 0.0_8
  ELSE
    rhxd = (qvd*qs-qv*qsd)/qs**2
    rhx = qv/qs
  END IF
  k1d = -(cons_alhl**2*rho_w*k_cond*cons_rvap*2*te*ted/(k_cond*cons_rvap&
&   *te**2)**2)
  k1 = cons_alhl**2*rho_w/(k_cond*cons_rvap*te**2)
  k2d = (cons_rvap*rho_w*ted*diffu*1000.*es/pl-cons_rvap*te*rho_w*diffu*&
&   1000.*esd/pl)/(diffu*(1000./pl)*es)**2
  k2 = cons_rvap*te*rho_w/(diffu*(1000./pl)*es)
!Here DIFFU is given for 1000 mb so 1000./PR accounts for increased diffusivity at lower pressure. 
  IF (f .GT. 0. .AND. ql .GT. 0.) THEN
    qcmd = (qld*f-ql*fd)/f**2
    qcm = ql/f
  ELSE
    qcm = 0.
    qcmd = 0.0_8
  END IF
qcmd = 0.0_8
  CALL ldradius_d(pl, te, ted, qcm, qcmd, nn, rho_w, radius, radiusd, &
&           cons_rgas, cons_pi)
  IF (rhx .LT. rhcr .AND. radius .GT. 0.0) THEN
! / (1.00 - RHx)
    teffd = (-(rhxd*(k1+k2)*radius**2)-(rhcr-rhx)*((k1d+k2d)*radius**2+(&
&     k1+k2)*2*radius*radiusd))/((k1+k2)*radius**2)**2
    teff = (rhcr-rhx)/((k1+k2)*radius**2)
  ELSE
! -999.
    teff = 0.0
    teffd = 0.0_8
  END IF
  evapd = a_eff*dt*(qld*teff+ql*teffd)
  evap = a_eff*ql*dt*teff
  IF (evap .GT. ql) THEN
    evapd = qld
    evap = ql
  ELSE
    evap = evap
  END IF
  qcd = qld + qid
  qc = ql + qi
  IF (qc .GT. 0.) THEN
    fd = ((fd*(qc-evap)+f*(qcd-evapd))*qc-f*(qc-evap)*qcd)/qc**2
    f = f*(qc-evap)/qc
  END IF
  qvd = qvd + evapd
  qv = qv + evap
  qld = qld - evapd
  ql = ql - evap
  ted = ted - cons_alhl*evapd/cons_cp
  te = te - cons_alhl/cons_cp*evap
END SUBROUTINE evap_cnv_d

!  Differentiation of subl_cnv in forward (tangent) mode:
!   variations   of useful results: f qi qv te
!   with respect to varying inputs: f qi ql qs qv te
SUBROUTINE subl_cnv_d(dt, rhcr, pl, te, ted, qv, qvd, ql, qld, qi, qid, &
& f, fd, xf, qs, qsd, rho_w, cld_evp_eff, cons_h2omw, cons_airmw, &
& cons_alhl, cons_rvap, cons_rgas, cons_pi, cons_cp, cons_alhs)
  IMPLICIT NONE
!INPUTS
  REAL*8, INTENT(IN) :: dt, rhcr, pl, xf, qs, rho_w, cld_evp_eff
  REAL*8, INTENT(IN) :: qsd
  REAL*8, INTENT(IN) :: cons_h2omw, cons_airmw, cons_alhl, cons_rvap, &
& cons_rgas, cons_pi, cons_cp, cons_alhs
!PROGNOSTIC
  REAL*8, INTENT(INOUT) :: te, qv, ql, qi, f
  REAL*8, INTENT(INOUT) :: ted, qvd, qld, qid, fd
!LOCALS
  REAL*8 :: es, radius, k1, k2, teff, qcm, subl, rhx, qc, a_eff, nn, &
& epsilon
  REAL*8 :: esd, radiusd, k1d, k2d, teffd, qcmd, subld, rhxd, qcd
  REAL*8, PARAMETER :: k_cond=2.4e-2
  REAL*8, PARAMETER :: diffu=2.2e-5
  INTRINSIC min
  epsilon = cons_h2omw/cons_airmw
  a_eff = cld_evp_eff
  nn = 5.*1.0e6
! (100 s <-^ convert from mbar to Pa)
  esd = (100.*pl*qsd*(epsilon+(1.0-epsilon)*qs)-100.*pl*qs*(1.0-epsilon)&
&   *qsd)/(epsilon+(1.0-epsilon)*qs)**2
  es = 100.*pl*qs/(epsilon+(1.0-epsilon)*qs)
  IF (qv/qs .GT. 1.00) THEN
    rhx = 1.00
    rhxd = 0.0_8
  ELSE
    rhxd = (qvd*qs-qv*qsd)/qs**2
    rhx = qv/qs
  END IF
  k1d = -(cons_alhl**2*rho_w*k_cond*cons_rvap*2*te*ted/(k_cond*cons_rvap&
&   *te**2)**2)
  k1 = cons_alhl**2*rho_w/(k_cond*cons_rvap*te**2)
  k2d = (cons_rvap*rho_w*ted*diffu*1000.*es/pl-cons_rvap*te*rho_w*diffu*&
&   1000.*esd/pl)/(diffu*(1000./pl)*es)**2
  k2 = cons_rvap*te*rho_w/(diffu*(1000./pl)*es)
!Here DIFFU is given for 1000 mb so 1000./PR accounts for increased diffusivity at lower pressure.
  IF (f .GT. 0. .AND. qi .GT. 0.) THEN
    qcmd = (qid*f-qi*fd)/f**2
    qcm = qi/f
  ELSE
    qcm = 0.
    qcmd = 0.0_8
  END IF
qcmd = 0.0_8
  CALL ldradius_d(pl, te, ted, qcm, qcmd, nn, rho_w, radius, radiusd, &
&           cons_rgas, cons_pi)
  IF (rhx .LT. rhcr .AND. radius .GT. 0.0) THEN
! / (1.00 - RHx)
    teffd = (-(rhxd*(k1+k2)*radius**2)-(rhcr-rhx)*((k1d+k2d)*radius**2+(&
&     k1+k2)*2*radius*radiusd))/((k1+k2)*radius**2)**2
    teff = (rhcr-rhx)/((k1+k2)*radius**2)
  ELSE
! -999.
    teff = 0.0
    teffd = 0.0_8
  END IF
  subld = a_eff*dt*(qid*teff+qi*teffd)
  subl = a_eff*qi*dt*teff
  IF (subl .GT. qi) THEN
    subld = qid
    subl = qi
  ELSE
    subl = subl
  END IF
  qcd = qld + qid
  qc = ql + qi
  IF (qc .GT. 0.) THEN
    fd = ((fd*(qc-subl)+f*(qcd-subld))*qc-f*(qc-subl)*qcd)/qc**2
    f = f*(qc-subl)/qc
  END IF
  qvd = qvd + subld
  qv = qv + subl
  qid = qid - subld
  qi = qi - subl
  ted = ted - cons_alhs*subld/cons_cp
  te = te - cons_alhs/cons_cp*subl
END SUBROUTINE subl_cnv_d

!  Differentiation of ldradius in forward (tangent) mode:
!   variations   of useful results: radius
!   with respect to varying inputs: qcl te
SUBROUTINE ldradius_d(pl, te, ted, qcl, qcld, nn, rho_w, radius, radiusd&
& , cons_rgas, cons_pi)
  IMPLICIT NONE
!Inputs
  REAL*8, INTENT(IN) :: te, pl, nn, qcl, rho_w
  REAL*8, INTENT(IN) :: ted, qcld
  REAL*8, INTENT(IN) :: cons_rgas, cons_pi
!Outputs      
  REAL*8, INTENT(OUT) :: radius
  REAL*8, INTENT(OUT) :: radiusd
  REAL*8 :: pwx1
  REAL*8 :: pwx1d
!Equiv. Spherical Cloud Particle Radius in m
  pwx1d = (qcld*100.*pl/(cons_rgas*te)-qcl*100.*pl*ted/(cons_rgas*te**2)&
&   )/(nn*rho_w*(4./3.)*cons_pi)
  pwx1 = qcl*(100.*pl/(cons_rgas*te))/(nn*rho_w*(4./3.)*cons_pi)
  IF (pwx1 .GT. 0.0 .OR. (pwx1 .LT. 0.0 .AND. 1./3. .EQ. int(1./3.))) &
& THEN
    radiusd = pwx1**(1./3.-1)*pwx1d/3.
  ELSE IF (pwx1 .EQ. 0.0 .AND. 1./3. .EQ. 1.0) THEN
    radiusd = pwx1d
  ELSE
    radiusd = 0.0
  END IF
  radius = pwx1**(1./3.)
END SUBROUTINE ldradius_d

!  Differentiation of autoconversion_ls in forward (tangent) mode:
!   variations   of useful results: f qc qp
!   with respect to varying inputs: f qc te
SUBROUTINE autoconversion_ls_d(dt, qc, qcd, qp, qpd, te, ted, pl, f, fd&
& , sundqv2, sundqv3, sundqt1, c_00, lwcrit, dzet)
  IMPLICIT NONE
!Inputs
  REAL*8, INTENT(IN) :: dt, te, pl, dzet, sundqv2, sundqv3, sundqt1, &
& c_00, lwcrit
  REAL*8, INTENT(IN) :: ted
!Prognostic
  REAL*8, INTENT(INOUT) :: qc, qp, f
  REAL*8, INTENT(INOUT) :: qcd, qpd, fd
!Locals
  REAL*8 :: acf0, acf, c00x, iqccrx, f2, f3, rate, dqp, qcm, dqfac
  REAL*8 :: c00xd, iqccrxd, f2d, f3d, rated, dqpd, qcmd, dqfacd
  INTRINSIC exp
  INTRINSIC min
  INTRINSIC max
  REAL*8 :: arg1
  REAL*8 :: arg1d
  REAL*8 :: x4
  REAL*8 :: x3
  REAL*8 :: x2
  REAL*8 :: x2d
  REAL*8 :: x1
  REAL*8 :: x1d
  REAL*8 :: x4d
!Zero Locals
  acf0 = 0.0
  acf = 0.0
  c00x = 0.0
  iqccrx = 0.0
  f2 = 0.0
  f3 = 0.0
  rate = 0.0
  dqp = 0.0
  qcm = 0.0
  dqfac = 0.0
  CALL cons_sundq3_d(te, ted, sundqv2, sundqv3, sundqt1, f2, f2d, f3)
  f2d = 0.5*f2d
  c00xd = c_00*f3*f2d
  c00x = c_00*f2*f3
  iqccrxd = f3*f2d/lwcrit
  iqccrx = f2*f3/lwcrit
  IF (f .GT. 0. .AND. qc .GT. 0.) THEN
    qcmd = (qcd*f-qc*fd)/f**2
    qcm = qc/f
  ELSE
    qcm = 0.
    qcmd = 0.0_8
  END IF
  arg1d = -(2*qcm*iqccrx*(qcmd*iqccrx+qcm*iqccrxd))
  arg1 = -((qcm*iqccrx)**2)
  rated = c00xd*(1.0-exp(arg1)) - c00x*arg1d*exp(arg1)
  rate = c00x*(1.0-exp(arg1))
!Temporary kluge until we can figure a better to make thicker low clouds.
  f2 = 1.0
  f3 = 1.0
!Implement ramps for gradual change in autoconv
!Thicken low high lat clouds
  IF (pl .GE. 775. .AND. te .LE. 275.) f3 = 0.2
!F3 = max(-0.016 * PL + 13.4, 0.2)
  IF (pl .GE. 825. .AND. te .LE. 282.) f3 = 0.2
!F3 = max(0.11 * TE - 30.02, 0.2)
  IF (pl .GE. 775. .AND. pl .LT. 825. .AND. te .LE. 282. .AND. te .GT. &
&     275.) f3 = 0.2
!F3 = min(max(-0.016*PL + 0.11 * TE - 16.85, 0.2),1.)
  IF (pl .GE. 825. .AND. te .LE. 275.) f3 = 0.2
  IF (pl .LE. 775. .OR. te .GT. 282.) f3 = 1.
!Thin-out low tropical clouds
  IF (pl .GE. 950. .AND. te .GE. 285.) THEN
    IF (0.2*te - 56 .GT. 2.) THEN
      f3 = 2.
      f3d = 0.0_8
    ELSE
      f3d = 0.2*ted
      f3 = 0.2*te - 56
    END IF
  ELSE
    f3d = 0.0_8
  END IF
  IF (pl .GE. 925. .AND. te .GE. 290.) THEN
    IF (0.04*pl - 36. .GT. 2.) THEN
      f3 = 2.
      f3d = 0.0_8
    ELSE
      f3 = 0.04*pl - 36.
      f3d = 0.0_8
    END IF
  END IF
  IF (pl .GE. 925. .AND. pl .LT. 950. .AND. te .GT. 285. .AND. te .LT. &
&     290.) THEN
    IF (0.04*pl + 0.2*te - 94. .GT. 2.) THEN
      x1 = 2.
      x1d = 0.0_8
    ELSE
      x1d = 0.2*ted
      x1 = 0.04*pl + 0.2*te - 94.
    END IF
    IF (x1 .LT. 1.) THEN
      f3 = 1.
      f3d = 0.0_8
    ELSE
      f3d = x1d
      f3 = x1
    END IF
  END IF
  IF (pl .GE. 950. .AND. te .GE. 290.) THEN
    f3 = 2.
    f3d = 0.0_8
  END IF
  IF (f3 .LT. 0.1) THEN
    f3 = 0.1
    f3d = 0.0_8
  ELSE
    f3 = f3
  END IF
  rated = f3d*rate + f3*rated
  rate = f3*rate
  dqpd = qcd*(1.0-exp(-(rate*dt))) + qc*dt*rated*exp(-(rate*dt))
  dqp = qc*(1.0-exp(-(rate*dt)))
  IF (dqp .LT. 0.0) THEN
    dqp = 0.0
    dqpd = 0.0_8
  ELSE
    dqp = dqp
  END IF
!Wipe-out warm fogs
  dqfac = 0.
  IF (pl .GE. 975. .AND. te .GE. 280.) THEN
    IF (0.2*te - 56. .GT. 1.) THEN
      x2 = 1.
      x2d = 0.0_8
    ELSE
      x2d = 0.2*ted
      x2 = 0.2*te - 56.
    END IF
    IF (x2 .LT. 0.) THEN
      dqfac = 0.
      dqfacd = 0.0_8
    ELSE
      dqfacd = x2d
      dqfac = x2
    END IF
  ELSE
    dqfacd = 0.0_8
  END IF
  IF (pl .GE. 950. .AND. te .GE. 285.) THEN
    IF (0.04*pl - 38. .GT. 1.) THEN
      x3 = 1.
    ELSE
      x3 = 0.04*pl - 38.
    END IF
    IF (x3 .LT. 0.) THEN
      dqfac = 0.
      dqfacd = 0.0_8
    ELSE
      dqfac = x3
      dqfacd = 0.0_8
    END IF
  END IF
  IF (pl .GE. 950. .AND. pl .LT. 975. .AND. te .GT. 280. .AND. te .LT. &
&     285.) THEN
    IF (0.04*pl + 0.2*te - 95. .GT. 1.) THEN
      x4 = 1.
      x4d = 0.0_8
    ELSE
      x4d = 0.2*ted
      x4 = 0.04*pl + 0.2*te - 95.
    END IF
    IF (x4 .LT. 0.) THEN
      dqfac = 0.
      dqfacd = 0.0_8
    ELSE
      dqfacd = x4d
      dqfac = x4
    END IF
  END IF
  IF (pl .GE. 975. .AND. te .GE. 285.) THEN
    dqfac = 1.
    dqfacd = 0.0_8
  END IF
  IF (dqp .LT. dqfac*qc) THEN
    dqpd = dqfacd*qc + dqfac*qcd
    dqp = dqfac*qc
  ELSE
    dqp = dqp
  END IF
  qcd = qcd - dqpd
  qc = qc - dqp
  qpd = dqpd
  qp = qp + dqp
!IF LARGE SCALE THEN
  IF (qc + dqp .GT. 0.) THEN
    fd = ((qcd*f+qc*fd)*(qc+dqp)-qc*f*(qcd+dqpd))/(qc+dqp)**2
    f = qc*f/(qc+dqp)
  END IF
END SUBROUTINE autoconversion_ls_d

!  Differentiation of autoconversion_cnv in forward (tangent) mode:
!   variations   of useful results: qc qp
!   with respect to varying inputs: f qc te
SUBROUTINE autoconversion_cnv_d(dt, qc, qcd, qp, qpd, te, ted, pl, f, fd&
& , sundqv2, sundqv3, sundqt1, c_00, lwcrit, dzet)
  IMPLICIT NONE
!Inputs
  REAL*8, INTENT(IN) :: dt, te, pl, dzet, sundqv2, sundqv3, sundqt1, &
& c_00, lwcrit
  REAL*8, INTENT(IN) :: ted
!Prognostic
  REAL*8, INTENT(INOUT) :: qc, qp, f
  REAL*8, INTENT(INOUT) :: qcd, qpd, fd
!Locals
  REAL*8 :: acf0, acf, c00x, iqccrx, f2, f3, rate, dqp, qcm, dqfac
  REAL*8 :: c00xd, iqccrxd, f2d, f3d, rated, dqpd, qcmd, dqfacd
  INTRINSIC exp
  INTRINSIC min
  INTRINSIC max
  REAL*8 :: arg1
  REAL*8 :: arg1d
  REAL*8 :: x4
  REAL*8 :: x3
  REAL*8 :: x2
  REAL*8 :: x2d
  REAL*8 :: x1
  REAL*8 :: x1d
  REAL*8 :: x4d
!Zero Locals
  acf0 = 0.0
  acf = 0.0
  c00x = 0.0
  iqccrx = 0.0
  f2 = 0.0
  f3 = 0.0
  rate = 0.0
  dqp = 0.0
  qcm = 0.0
  dqfac = 0.0
  CALL cons_sundq3_d(te, ted, sundqv2, sundqv3, sundqt1, f2, f2d, f3)
  f2d = 0.5*f2d
  c00xd = c_00*f3*f2d
  c00x = c_00*f2*f3
  iqccrxd = f3*f2d/lwcrit
  iqccrx = f2*f3/lwcrit
  IF (f .GT. 0. .AND. qc .GT. 0.) THEN
    qcmd = (qcd*f-qc*fd)/f**2
    qcm = qc/f
  ELSE
    qcm = 0.
    qcmd = 0.0_8
  END IF
  arg1d = -(2*qcm*iqccrx*(qcmd*iqccrx+qcm*iqccrxd))
  arg1 = -((qcm*iqccrx)**2)
  rated = c00xd*(1.0-exp(arg1)) - c00x*arg1d*exp(arg1)
  rate = c00x*(1.0-exp(arg1))
!Temporary kluge until we can figure a better to make thicker low clouds.
  f2 = 1.0
  f3 = 1.0
!Implement ramps for gradual change in autoconv
!Thicken low high lat clouds
  IF (pl .GE. 775. .AND. te .LE. 275.) f3 = 0.2
!F3 = max(-0.016 * PL + 13.4, 0.2)
  IF (pl .GE. 825. .AND. te .LE. 282.) f3 = 0.2
!F3 = max(0.11 * TE - 30.02, 0.2)
  IF (pl .GE. 775. .AND. pl .LT. 825. .AND. te .LE. 282. .AND. te .GT. &
&     275.) f3 = 0.2
!F3 = min(max(-0.016*PL + 0.11 * TE - 16.85, 0.2),1.)
  IF (pl .GE. 825. .AND. te .LE. 275.) f3 = 0.2
  IF (pl .LE. 775. .OR. te .GT. 282.) f3 = 1.
!Thin-out low tropical clouds
  IF (pl .GE. 950. .AND. te .GE. 285.) THEN
    IF (0.2*te - 56 .GT. 2.) THEN
      f3 = 2.
      f3d = 0.0_8
    ELSE
      f3d = 0.2*ted
      f3 = 0.2*te - 56
    END IF
  ELSE
    f3d = 0.0_8
  END IF
  IF (pl .GE. 925. .AND. te .GE. 290.) THEN
    IF (0.04*pl - 36. .GT. 2.) THEN
      f3 = 2.
      f3d = 0.0_8
    ELSE
      f3 = 0.04*pl - 36.
      f3d = 0.0_8
    END IF
  END IF
  IF (pl .GE. 925. .AND. pl .LT. 950. .AND. te .GT. 285. .AND. te .LT. &
&     290.) THEN
    IF (0.04*pl + 0.2*te - 94. .GT. 2.) THEN
      x1 = 2.
      x1d = 0.0_8
    ELSE
      x1d = 0.2*ted
      x1 = 0.04*pl + 0.2*te - 94.
    END IF
    IF (x1 .LT. 1.) THEN
      f3 = 1.
      f3d = 0.0_8
    ELSE
      f3d = x1d
      f3 = x1
    END IF
  END IF
  IF (pl .GE. 950. .AND. te .GE. 290.) THEN
    f3 = 2.
    f3d = 0.0_8
  END IF
  IF (f3 .LT. 0.1) THEN
    f3 = 0.1
    f3d = 0.0_8
  ELSE
    f3 = f3
  END IF
  rated = f3d*rate + f3*rated
  rate = f3*rate
  dqpd = qcd*(1.0-exp(-(rate*dt))) + qc*dt*rated*exp(-(rate*dt))
  dqp = qc*(1.0-exp(-(rate*dt)))
  IF (dqp .LT. 0.0) THEN
    dqp = 0.0
    dqpd = 0.0_8
  ELSE
    dqp = dqp
  END IF
!Wipe-out warm fogs
  dqfac = 0.
  IF (pl .GE. 975. .AND. te .GE. 280.) THEN
    IF (0.2*te - 56. .GT. 1.) THEN
      x2 = 1.
      x2d = 0.0_8
    ELSE
      x2d = 0.2*ted
      x2 = 0.2*te - 56.
    END IF
    IF (x2 .LT. 0.) THEN
      dqfac = 0.
      dqfacd = 0.0_8
    ELSE
      dqfacd = x2d
      dqfac = x2
    END IF
  ELSE
    dqfacd = 0.0_8
  END IF
  IF (pl .GE. 950. .AND. te .GE. 285.) THEN
    IF (0.04*pl - 38. .GT. 1.) THEN
      x3 = 1.
    ELSE
      x3 = 0.04*pl - 38.
    END IF
    IF (x3 .LT. 0.) THEN
      dqfac = 0.
      dqfacd = 0.0_8
    ELSE
      dqfac = x3
      dqfacd = 0.0_8
    END IF
  END IF
  IF (pl .GE. 950. .AND. pl .LT. 975. .AND. te .GT. 280. .AND. te .LT. &
&     285.) THEN
    IF (0.04*pl + 0.2*te - 95. .GT. 1.) THEN
      x4 = 1.
      x4d = 0.0_8
    ELSE
      x4d = 0.2*ted
      x4 = 0.04*pl + 0.2*te - 95.
    END IF
    IF (x4 .LT. 0.) THEN
      dqfac = 0.
      dqfacd = 0.0_8
    ELSE
      dqfacd = x4d
      dqfac = x4
    END IF
  END IF
  IF (pl .GE. 975. .AND. te .GE. 285.) THEN
    dqfac = 1.
    dqfacd = 0.0_8
  END IF
  IF (dqp .LT. dqfac*qc) THEN
    dqpd = dqfacd*qc + dqfac*qcd
    dqp = dqfac*qc
  ELSE
    dqp = dqp
  END IF
  qcd = qcd - dqpd
  qc = qc - dqp
  qpd = dqpd
  qp = qp + dqp
END SUBROUTINE autoconversion_cnv_d

!  Differentiation of get_ice_fraction in forward (tangent) mode:
!   variations   of useful results: icefrct
!   with respect to varying inputs: temp
SUBROUTINE get_ice_fraction_d(temp, tempd, t_ice_all, t_ice_max, &
& icefrpwr, icefrct, icefrctd)
  IMPLICIT NONE
!Inputs
  REAL*8, INTENT(IN) :: temp, t_ice_all, t_ice_max
  REAL*8, INTENT(IN) :: tempd
  INTEGER, INTENT(IN) :: icefrpwr
!Outputs
  REAL*8, INTENT(OUT) :: icefrct
  REAL*8, INTENT(OUT) :: icefrctd
  INTRINSIC min
  INTRINSIC max
  icefrct = 0.00
  IF (temp .LE. t_ice_all) THEN
    icefrct = 1.000
    icefrctd = 0.0_8
  ELSE IF (temp .GT. t_ice_all .AND. temp .LE. t_ice_max) THEN
    icefrctd = -(tempd/(t_ice_max-t_ice_all))
    icefrct = 1.00 - (temp-t_ice_all)/(t_ice_max-t_ice_all)
  ELSE
    icefrctd = 0.0_8
  END IF
  IF (icefrct .GT. 1.00) THEN
    icefrct = 1.00
    icefrctd = 0.0_8
  ELSE
    icefrct = icefrct
  END IF
  IF (icefrct .LT. 0.00) THEN
    icefrct = 0.00
    icefrctd = 0.0_8
  ELSE
    icefrct = icefrct
  END IF
  IF (icefrct .GT. 0.0 .OR. (icefrct .LT. 0.0 .AND. icefrpwr .EQ. int(&
&     icefrpwr))) THEN
    icefrctd = icefrpwr*icefrct**(icefrpwr-1)*icefrctd
  ELSE IF (icefrct .EQ. 0.0 .AND. icefrpwr .EQ. 1.0) THEN
    icefrctd = icefrctd
  ELSE
    icefrctd = 0.0
  END IF
  icefrct = icefrct**icefrpwr
END SUBROUTINE get_ice_fraction_d

!  Differentiation of cons_sundq3 in forward (tangent) mode:
!   variations   of useful results: f2
!   with respect to varying inputs: temp
SUBROUTINE cons_sundq3_d(temp, tempd, rate2, rate3, te1, f2, f2d, f3)
  IMPLICIT NONE
!Inputs
  REAL*8, INTENT(IN) :: rate2, rate3, te1, temp
  REAL*8, INTENT(IN) :: tempd
!Outputs
  REAL*8, INTENT(OUT) :: f2, f3
  REAL*8, INTENT(OUT) :: f2d
!Locals
  REAL*8, PARAMETER :: te0=273.
  REAL*8, PARAMETER :: te2=200.
  REAL*8 :: jump1
  INTRINSIC abs
  INTRINSIC min
  REAL*8 :: abs0
  jump1 = (rate2-1.0)/(te0-te1)**0.333
!Ice - phase treatment 
  IF (temp .GE. te0) THEN
    f2 = 1.0
    f3 = 1.0
  END IF
  IF (temp .GE. te1 .AND. temp .LT. te0) THEN
    IF (te0 - temp .GE. 0.) THEN
      abs0 = te0 - temp
    ELSE
      abs0 = -(te0-temp)
    END IF
    IF (abs0 .GT. 0.0) THEN
!Linearisation security
      f2d = -(jump1*0.3333*(te0-temp)**(-0.6667)*tempd)
      f2 = 1.0 + jump1*(te0-temp)**0.3333
    ELSE
      f2 = 1.0
      f2d = 0.0_8
    END IF
    f3 = 1.0
  ELSE
    f2d = 0.0_8
  END IF
  IF (temp .LT. te1) THEN
    f2d = (-((rate3-rate2)*tempd))/(te1-te2)
    f2 = rate2 + (rate3-rate2)*(te1-temp)/(te1-te2)
    f3 = 1.0
  END IF
  IF (f2 .GT. 27.0) THEN
    f2 = 27.0
    f2d = 0.0_8
  ELSE
    f2 = f2
  END IF
END SUBROUTINE cons_sundq3_d

!  Differentiation of cons_microphys in forward (tangent) mode:
!   variations   of useful results: aa bb
!   with respect to varying inputs: temp q_sat alhx3
SUBROUTINE cons_microphys_d(temp, tempd, pr, q_sat, q_satd, aa, aad, bb&
& , bbd, cons_h2omw, cons_airmw, cons_rvap, alhx3, alhx3d)
  IMPLICIT NONE
!Inputs
  REAL*8, INTENT(IN) :: temp, q_sat, pr, alhx3
  REAL*8, INTENT(IN) :: tempd, q_satd, alhx3d
  REAL*8, INTENT(IN) :: cons_h2omw, cons_airmw, cons_rvap
!Outputs
  REAL*8, INTENT(OUT) :: aa, bb
  REAL*8, INTENT(OUT) :: aad, bbd
!Locals
  REAL*8, PARAMETER :: k_cond=2.4e-2
  REAL*8, PARAMETER :: diffu=2.2e-5
  REAL*8 :: e_sat, epsi
  REAL*8 :: e_satd
  epsi = cons_h2omw/cons_airmw
! (100 converts from mbar to Pa)
  e_satd = (100.*pr*q_satd*(epsi+(1.0-epsi)*q_sat)-100.*pr*q_sat*(1.0-&
&   epsi)*q_satd)/(epsi+(1.0-epsi)*q_sat)**2
  e_sat = 100.*pr*q_sat/(epsi+(1.0-epsi)*q_sat)
  aad = (2*alhx3*alhx3d*k_cond*cons_rvap*temp**2-alhx3**2*k_cond*&
&   cons_rvap*2*temp*tempd)/(k_cond*cons_rvap*temp**2)**2
  aa = alhx3**2/(k_cond*cons_rvap*temp**2)
  bbd = (cons_rvap*tempd*diffu*1000.*e_sat/pr-cons_rvap*temp*diffu*1000.&
&   *e_satd/pr)/(diffu*(1000./pr)*e_sat)**2
  bb = cons_rvap*temp/(diffu*(1000./pr)*e_sat)
END SUBROUTINE cons_microphys_d

!  Differentiation of cons_alhx in forward (tangent) mode:
!   variations   of useful results: alhx3
!   with respect to varying inputs: t alhx3
SUBROUTINE cons_alhx_d(t, td, alhx3, alhx3d, t_ice_max, t_ice_all, &
& cons_alhs, cons_alhl)
  IMPLICIT NONE
!Inputs
  REAL*8, INTENT(IN) :: t, t_ice_max, t_ice_all
  REAL*8, INTENT(IN) :: td
  REAL*8, INTENT(IN) :: cons_alhs, cons_alhl
!Outputs
  REAL*8, INTENT(OUT) :: alhx3
  REAL*8, INTENT(OUT) :: alhx3d
  IF (t .LT. t_ice_all) THEN
    alhx3 = cons_alhs
    alhx3d = 0.0_8
  END IF
  IF (t .GT. t_ice_max) THEN
    alhx3 = cons_alhl
    alhx3d = 0.0_8
  END IF
  IF (t .LE. t_ice_max .AND. t .GE. t_ice_all) THEN
    alhx3d = (cons_alhl-cons_alhs)*td/(t_ice_max-t_ice_all)
    alhx3 = cons_alhs + (cons_alhl-cons_alhs)*(t-t_ice_all)/(t_ice_max-&
&     t_ice_all)
  END IF
END SUBROUTINE cons_alhx_d

!  Differentiation of ice_settlefall_cnv in forward (tangent) mode:
!   variations   of useful results: qi qp
!   with respect to varying inputs: f qi dz te
SUBROUTINE ice_settlefall_cnv_d(wxr, qi, qid, pl, te, ted, f, fd, &
& cons_rgas, khu, khl, k, dt, dz, dzd, qp, qpd, anv_icefall_c)
  IMPLICIT NONE
!Inputs
  REAL*8, INTENT(IN) :: wxr, pl, te, dz, dt, anv_icefall_c
  REAL*8, INTENT(IN) :: ted, dzd
  REAL*8, INTENT(IN) :: cons_rgas
  INTEGER, INTENT(IN) :: khu, khl, k
  REAL*8, INTENT(INOUT) :: qi, f, qp
  REAL*8, INTENT(INOUT) :: qid, fd, qpd
!Locals
  REAL*8 :: rho, xim, lxim, qixp, vf
  REAL*8 :: rhod, ximd, lximd, qixpd, vfd
  INTRINSIC log10
  INTRINSIC max
  INTRINSIC min
  REAL*8 :: pwx1
  REAL*8 :: pwr1
  REAL*8 :: max1
! 1000 TAKES TO g m^-3 ; 100 takes mb TO Pa
  rhod = -(1000.*100.*pl*cons_rgas*ted/(cons_rgas*te)**2)
  rho = 1000.*100.*pl/(cons_rgas*te)
  IF (f .GT. 0. .AND. qi .GT. 0.) THEN
    ximd = (qid*f-qi*fd)*rho/f**2 + qi*rhod/f
    xim = qi/f*rho
  ELSE
    xim = 0.
    ximd = 0.0_8
  END IF
  IF (xim .GT. 0.) THEN
    lximd = ximd/(xim*log(10.0))
    lxim = log10(xim)
  ELSE
    lxim = 0.0
    lximd = 0.0_8
  END IF
  vfd = 53.2*lximd + 5.5*2*lxim*lximd
  vf = 128.6 + 53.2*lxim + 5.5*lxim**2
!VF = VF*100./MAX(PL,10.) ! Reduce/increase fall speeds for high/low pressure (NOT in LC98!!! ) 
! Assume unmodified they represent situation at 100 mb
  IF (wxr .GT. 0.) THEN
    IF (pl .LT. 10.) THEN
      max1 = 10.
    ELSE
      max1 = pl
    END IF
    pwx1 = 100./max1
    pwr1 = pwx1**wxr
    vfd = pwr1*vfd
    vf = vf*pwr1
  END IF
  vfd = vfd/100.
  vf = vf/100.
  IF (khu .GT. 0 .AND. khl .GT. 0) THEN
    IF (k - 1 .GE. khu .AND. k - 1 .LE. khl) THEN
      vfd = 0.01*vfd
      vf = 0.01*vf
    END IF
  END IF
  vfd = anv_icefall_c*vfd
  vf = anv_icefall_c*vf
  qixp = 0.0
  qixpd = qid*vf*dt/dz + qi*(dt*vfd*dz-vf*dt*dzd)/dz**2
  qixp = qi*(vf*dt/dz)
  IF (qixp .GT. qi) THEN
    qixpd = qid
    qixp = qi
  ELSE
    qixp = qixp
  END IF
  IF (qixp .LT. 0.0) THEN
    qixp = 0.0
    qixpd = 0.0_8
  ELSE
    qixp = qixp
  END IF
  qpd = qixpd
  qp = qp + qixp
  qid = qid - qixpd
  qi = qi - qixp
END SUBROUTINE ice_settlefall_cnv_d

!  Differentiation of ice_settlefall_ls in forward (tangent) mode:
!   variations   of useful results: f qi qp
!   with respect to varying inputs: f qi dz te
SUBROUTINE ice_settlefall_ls_d(wxr, qi, qid, pl, te, ted, f, fd, &
& cons_rgas, khu, khl, k, dt, dz, dzd, qp, qpd, ls_icefall_c)
  IMPLICIT NONE
!Inputs
  REAL*8, INTENT(IN) :: wxr, pl, te, dz, dt, ls_icefall_c
  REAL*8, INTENT(IN) :: ted, dzd
  REAL*8, INTENT(IN) :: cons_rgas
  INTEGER, INTENT(IN) :: khu, khl, k
  REAL*8, INTENT(INOUT) :: qi, f, qp
  REAL*8, INTENT(INOUT) :: qid, fd, qpd
!Locals
  REAL*8 :: rho, xim, lxim, qixp, vf
  REAL*8 :: rhod, ximd, qixpd, vfd
  INTRINSIC log10
  INTRINSIC abs
  INTRINSIC max
  INTRINSIC min
  REAL*8 :: pwx1
  REAL*8 :: pwr1
  REAL*8 :: abs0
  REAL*8 :: max1
! 1000 TAKES TO g m^-3 ; 100 takes mb TO Pa
  rhod = -(1000.*100.*pl*cons_rgas*ted/(cons_rgas*te)**2)
  rho = 1000.*100.*pl/(cons_rgas*te)
  IF (f .GT. 0. .AND. qi .GT. 0.) THEN
    ximd = (qid*f-qi*fd)*rho/f**2 + qi*rhod/f
    xim = qi/f*rho
  ELSE
    xim = 0.
    ximd = 0.0_8
  END IF
  IF (xim .GT. 0.) THEN
    lxim = log10(xim)
  ELSE
    lxim = 0.0
  END IF
  IF (xim .GE. 0.) THEN
    abs0 = xim
  ELSE
    abs0 = -xim
  END IF
  IF (abs0 .GT. 0.0) THEN
!Linearisation security
    vfd = 109.0*0.16*xim**(-0.84)*ximd
    vf = 109.0*xim**0.16
  ELSE
    vf = 0.0
    vfd = 0.0_8
  END IF
!VF = VF*100./MAX(PL,10.) ! Reduce/increase fall speeds for high/low pressure (NOT in LC98!!! ) 
! Assume unmodified they represent situation at 100 mb
  IF (wxr .GT. 0.) THEN
    IF (pl .LT. 10.) THEN
      max1 = 10.
    ELSE
      max1 = pl
    END IF
    pwx1 = 100./max1
    pwr1 = pwx1**wxr
    vfd = pwr1*vfd
    vf = vf*pwr1
  END IF
  vfd = vfd/100.
  vf = vf/100.
  IF (khu .GT. 0 .AND. khl .GT. 0) THEN
    IF (k - 1 .GE. khu .AND. k - 1 .LE. khl) THEN
      vfd = 0.01*vfd
      vf = 0.01*vf
    END IF
  END IF
  vfd = ls_icefall_c*vfd
  vf = ls_icefall_c*vf
  qixp = 0.0
  qixpd = qid*vf*dt/dz + qi*(dt*vfd*dz-vf*dt*dzd)/dz**2
  qixp = qi*(vf*dt/dz)
  IF (qixp .GT. qi) THEN
    qixpd = qid
    qixp = qi
  ELSE
    qixp = qixp
  END IF
  IF (qixp .LT. 0.0) THEN
    qixp = 0.0
    qixpd = 0.0_8
  ELSE
    qixp = qixp
  END IF
  qpd = qixpd
  qp = qp + qixp
  qid = qid - qixpd
  qi = qi - qixp
  IF (qi + qixp .GT. 0.) THEN
    fd = ((qid*f+qi*fd)*(qi+qixp)-qi*f*(qid+qixpd))/(qi+qixp)**2
    f = qi*f/(qi+qixp)
  END IF
END SUBROUTINE ice_settlefall_ls_d

!  Differentiation of precipandevap in forward (tangent) mode:
!   variations   of useful results: evap_dd_above_out qv subl_dd_above_out
!                qcl pfi_above_out pfl_above_out te
!   with respect to varying inputs: aa area pfl_above_in qv pfi_above_in
!                bb evap_dd_above_in qcl qpi qpl subl_dd_above_in
!                dze qddf3 te
SUBROUTINE precipandevap_d(k, ktop, lm, dt, frland, rhcr3, qpl, qpld, &
& qpi, qpid, qcl, qcld, qci, te, ted, qv, qvd, mass, imass, pl, dze, &
& dzed, qddf3, qddf3d, aa, aad, bb, bbd, area, aread, pfl_above_in, &
& pfl_above_ind, pfl_above_out, pfl_above_outd, pfi_above_in, &
& pfi_above_ind, pfi_above_out, pfi_above_outd, evap_dd_above_in, &
& evap_dd_above_ind, evap_dd_above_out, evap_dd_above_outd, &
& subl_dd_above_in, subl_dd_above_ind, subl_dd_above_out, &
& subl_dd_above_outd, envfc, ddrfc, cons_alhf, cons_alhs, cons_alhl, &
& cons_cp, cons_tice, cons_h2omw, cons_airmw, revap_off_p, c_acc, c_ev_r&
& , c_ev_s, rho_w, estblx)
  IMPLICIT NONE
!Inputs
  INTEGER, INTENT(IN) :: k, lm, ktop
  REAL*8, INTENT(IN) :: dt, mass, imass, pl, aa, bb, rhcr3, dze, qddf3, &
& area, frland, envfc, ddrfc
  REAL*8, INTENT(IN) :: aad, bbd, dzed, qddf3d, aread
  REAL*8, INTENT(IN) :: cons_alhf, cons_alhs, cons_alhl, cons_cp, &
& cons_tice, cons_h2omw, cons_airmw
  REAL*8, INTENT(IN) :: revap_off_p
  REAL*8, INTENT(IN) :: c_acc, c_ev_r, c_ev_s, rho_w
  REAL*8, INTENT(IN) :: estblx(:)
!Prognostics
  REAL*8, INTENT(INOUT) :: qv, qpl, qpi, qcl, qci, te
  REAL*8, INTENT(INOUT) :: qvd, qpld, qpid, qcld, ted
  REAL*8, INTENT(INOUT) :: pfl_above_in, pfl_above_out, pfi_above_in, &
& pfi_above_out
  REAL*8, INTENT(INOUT) :: pfl_above_ind, pfl_above_outd, pfi_above_ind&
& , pfi_above_outd
  REAL*8, INTENT(INOUT) :: evap_dd_above_in, evap_dd_above_out, &
& subl_dd_above_in, subl_dd_above_out
  REAL*8, INTENT(INOUT) :: evap_dd_above_ind, evap_dd_above_outd, &
& subl_dd_above_ind, subl_dd_above_outd
!Locals
  INTEGER :: ns, nsmx, itr, l
  REAL*8 :: pfi, pfl, qs, dqs, envfrac, tko, qko, qstko, dqstko, rh_box&
& , t_ed, qplko, qpiko
  REAL*8 :: pfid, pfld, qsd, dqsd, tkod, qkod, qstkod, dqstkod, rh_boxd&
& , t_edd
  REAL*8 :: ifactor, rainrat0, snowrat0, fallrn, fallsn, vesn, vern, &
& nrain, nsnow, efactor
  REAL*8 :: ifactord, rainrat0d, snowrat0d, fallrnd, fallsnd, vesnd, &
& vernd, efactord
  REAL*8 :: tinlayerrn, diamrn, droprad, tinlayersn, diamsn, flakrad
  REAL*8 :: tinlayerrnd, diamrnd, dropradd, tinlayersnd, diamsnd, &
& flakradd
  REAL*8 :: evap, subl, accr, mltfrz, evapx, sublx, evap_dd, subl_dd, &
& ddfract, landseaf
  REAL*8 :: evapd, subld, accrd, mltfrzd, evap_ddd, subl_ddd
  REAL*8 :: tau_frz, tau_mlt
!m/s
  REAL*8, PARAMETER :: trmv_l=1.0
  LOGICAL, PARAMETER :: taneff=.false.
!Fraction of precip falling through "environment" vs through cloud
  REAL*8, PARAMETER :: b_sub=1.00
  INTRINSIC max
  INTRINSIC min
  INTRINSIC exp
  REAL*8 :: arg1
  REAL*8 :: arg1d
  envfrac = envfc
  IF (area .GT. 0.) THEN
    ifactord = -(aread/area**2)
    ifactor = 1./area
  ELSE
    ifactor = 1.00
    ifactord = 0.0_8
  END IF
  IF (ifactor .LT. 1.) THEN
    ifactor = 1.
    ifactord = 0.0_8
  ELSE
    ifactor = ifactor
  END IF
!Start at top of precip column:
!
!   a) Accrete                   
!   b) Evaporate/Sublimate  
!   c) Rain/Snow-out to next level down 
!   d) return to (a)
!Update saturated humidity
  CALL dqsats_bac_d(dqs, dqsd, qs, qsd, te, ted, pl, estblx, cons_h2omw&
&             , cons_airmw)
  ddfract = ddrfc
  IF (k .EQ. ktop) THEN
    pfld = mass*qpld
    pfl = qpl*mass
    pfid = mass*qpid
    pfi = qpi*mass
    evap_dd = 0.
    subl_dd = 0.
    evap_ddd = 0.0_8
    subl_ddd = 0.0_8
  ELSE
    qpld = qpld + imass*pfl_above_ind
    qpl = qpl + pfl_above_in*imass
    pfl = 0.00
    qpid = qpid + imass*pfi_above_ind
    qpi = qpi + pfi_above_in*imass
    pfi = 0.00
    accrd = b_sub*c_acc*mass*(qpld*qcl+qpl*qcld)
    accr = b_sub*c_acc*(qpl*mass)*qcl
    IF (accr .GT. qcl) THEN
      accrd = qcld
      accr = qcl
    ELSE
      accr = accr
    END IF
    qpld = qpld + accrd
    qpl = qpl + accr
    qcld = qcld - accrd
    qcl = qcl - accr
!Accretion of liquid condensate by falling ice/snow
    accrd = b_sub*c_acc*mass*(qpid*qcl+qpi*qcld)
    accr = b_sub*c_acc*(qpi*mass)*qcl
    IF (accr .GT. qcl) THEN
      accrd = qcld
      accr = qcl
    ELSE
      accr = accr
    END IF
    qpid = qpid + accrd
    qpi = qpi + accr
    qcld = qcld - accrd
    qcl = qcl - accr
!! Liquid freezes when accreted by snow
    ted = ted + cons_alhf*accrd/cons_cp
    te = te + cons_alhf*accr/cons_cp
    rainrat0d = mass*(ifactord*qpl+ifactor*qpld)/dt
    rainrat0 = ifactor*qpl*mass/dt
    snowrat0d = mass*(ifactord*qpi+ifactor*qpid)/dt
    snowrat0 = ifactor*qpi*mass/dt
    CALL marshpalm_d(rainrat0, rainrat0d, pl, diamrn, diamrnd, nrain, &
&              fallrn, fallrnd, vern, vernd)
    CALL marshpalm_d(snowrat0, snowrat0d, pl, diamsn, diamsnd, nsnow, &
&              fallsn, fallsnd, vesn, vesnd)
    tinlayerrnd = (dzed*(fallrn+0.01)-dze*fallrnd)/(fallrn+0.01)**2
    tinlayerrn = dze/(fallrn+0.01)
    tinlayersnd = (dzed*(fallsn+0.01)-dze*fallsnd)/(fallsn+0.01)**2
    tinlayersn = dze/(fallsn+0.01)
!Melting of Frozen precipitation      
! time scale for freezing (s). 
    tau_frz = 5000.
    mltfrz = 0.0
    IF (te .GT. cons_tice .AND. te .LE. cons_tice + 5.) THEN
      mltfrzd = ((tinlayersnd*qpi+tinlayersn*qpid)*(te-cons_tice)+&
&       tinlayersn*qpi*ted)/tau_frz
      mltfrz = tinlayersn*qpi*(te-cons_tice)/tau_frz
      IF (qpi .GT. mltfrz) THEN
        mltfrz = mltfrz
      ELSE
        mltfrzd = qpid
        mltfrz = qpi
      END IF
      ted = ted - cons_alhf*mltfrzd/cons_cp
      te = te - cons_alhf*mltfrz/cons_cp
      qpld = qpld + mltfrzd
      qpl = qpl + mltfrz
      qpid = qpid - mltfrzd
      qpi = qpi - mltfrz
    END IF
    mltfrz = 0.0
    IF (te .GT. cons_tice + 5.) THEN
! Go Ahead and melt any snow/hail left above 5 C 
      mltfrzd = qpid
      mltfrz = qpi
      ted = ted - cons_alhf*mltfrzd/cons_cp
      te = te - cons_alhf*mltfrz/cons_cp
      qpld = qpld + mltfrzd
      qpl = qpl + mltfrz
      qpid = qpid - mltfrzd
      qpi = qpi - mltfrz
    END IF
    mltfrz = 0.0
    IF (k .GE. lm - 1) THEN
      IF (te .GT. cons_tice + 0.) THEN
! Go Ahead and melt any snow/hail left above 0 C in lowest layers 
        mltfrzd = qpid
        mltfrz = qpi
        ted = ted - cons_alhf*mltfrzd/cons_cp
        te = te - cons_alhf*mltfrz/cons_cp
        qpld = qpld + mltfrzd
        qpl = qpl + mltfrz
        qpid = qpid - mltfrzd
        qpi = qpi - mltfrz
      END IF
    END IF
!Freezing of liquid precipitation      
    mltfrz = 0.0
    IF (te .LE. cons_tice) THEN
      ted = ted + cons_alhf*qpld/cons_cp
      te = te + cons_alhf*qpl/cons_cp
      qpid = qpld + qpid
      qpi = qpl + qpi
      mltfrz = qpl
      qpl = 0.
      qpld = 0.0_8
    END IF
!In the exp below, evaporation time scale is determined "microphysically" from temp, 
!press, and drop size. In this context C_EV becomes a dimensionless fudge-fraction. 
!Also remember that these microphysics are still only for liquid.
    qkod = qvd
    qko = qv
    tkod = ted
    tko = te
    qplko = qpl
    qpiko = qpi
!do itr = 1,1
    itr = 1
    dqstkod = dqsd
    dqstko = dqs
    qstkod = qsd + dqstkod*(tko-te) + dqstko*(tkod-ted)
    qstko = qs + dqstko*(tko-te)
    IF (qstko .LT. 1.0e-7) THEN
      qstko = 1.0e-7
      qstkod = 0.0_8
    ELSE
      qstko = qstko
    END IF
    rh_boxd = (qkod*qstko-qko*qstkod)/qstko**2
    rh_box = qko/qstko
    qko = qv
    tko = te
    IF (rh_box .LT. rhcr3) THEN
      efactord = (rho_w*(aad+bbd)*(rhcr3-rh_box)+rho_w*(aa+bb)*rh_boxd)/&
&       (rhcr3-rh_box)**2
      efactor = rho_w*(aa+bb)/(rhcr3-rh_box)
    ELSE
      efactor = 9.99e9
      efactord = 0.0_8
    END IF
    IF (frland .LT. 0.1) THEN
! Over Ocean
      landseaf = 0.5
    ELSE
! Over Land
      landseaf = 0.5
    END IF
    landseaf = 1.00
!Rain falling
    IF (rh_box .LT. rhcr3 .AND. diamrn .GT. 0.00 .AND. pl .GT. 100. &
&       .AND. pl .LT. revap_off_p) THEN
      dropradd = 0.5*diamrnd
      droprad = 0.5*diamrn
      t_edd = efactord*droprad**2 + efactor*2*droprad*dropradd
      t_ed = efactor*droprad**2
      t_edd = t_edd*(1.0+dqstko*cons_alhl/cons_cp) + t_ed*cons_alhl*&
&       dqstkod/cons_cp
      t_ed = t_ed*(1.0+dqstko*cons_alhl/cons_cp)
      arg1d = -((c_ev_r*landseaf*envfrac*(vernd*tinlayerrn+vern*&
&       tinlayerrnd)*t_ed-c_ev_r*vern*landseaf*envfrac*tinlayerrn*t_edd)&
&       /t_ed**2)
      arg1 = -(c_ev_r*vern*landseaf*envfrac*tinlayerrn/t_ed)
      evapd = qpld*(1.0-exp(arg1)) - qpl*arg1d*exp(arg1)
      evap = qpl*(1.0-exp(arg1))
    ELSE
      evap = 0.0
      evapd = 0.0_8
    END IF
!Snow falling
    IF (rh_box .LT. rhcr3 .AND. diamsn .GT. 0.00 .AND. pl .GT. 100. &
&       .AND. pl .LT. revap_off_p) THEN
      flakradd = 0.5*diamsnd
      flakrad = 0.5*diamsn
      t_edd = efactord*flakrad**2 + efactor*2*flakrad*flakradd
      t_ed = efactor*flakrad**2
      t_edd = t_edd*(1.0+dqstko*cons_alhs/cons_cp) + t_ed*cons_alhs*&
&       dqstkod/cons_cp
      t_ed = t_ed*(1.0+dqstko*cons_alhs/cons_cp)
      arg1d = -((c_ev_s*landseaf*envfrac*(vesnd*tinlayersn+vesn*&
&       tinlayersnd)*t_ed-c_ev_s*vesn*landseaf*envfrac*tinlayersn*t_edd)&
&       /t_ed**2)
      arg1 = -(c_ev_s*vesn*landseaf*envfrac*tinlayersn/t_ed)
      subld = qpid*(1.0-exp(arg1)) - qpi*arg1d*exp(arg1)
      subl = qpi*(1.0-exp(arg1))
    ELSE
      subl = 0.0
      subld = 0.0_8
    END IF
!if (itr == 1) then 
!   EVAPx  = EVAP
!   SUBLx  = SUBL
!else
!   EVAP   = (EVAP+EVAPx) /2.0
!   SUBL   = (SUBL+SUBLx) /2.0
!endif
    qko = qv + evap + subl
    tko = te - evap*cons_alhl/cons_cp - subl*cons_alhs/cons_cp
!enddo
    qpid = qpid - subld
    qpi = qpi - subl
    qpld = qpld - evapd
    qpl = qpl - evap
!Put some re-evap/re-subl precip in to a \quote{downdraft} to be applied later
    evap_ddd = evap_dd_above_ind + ddfract*mass*evapd
    evap_dd = evap_dd_above_in + ddfract*evap*mass
    evapd = evapd - ddfract*evapd
    evap = evap - ddfract*evap
    subl_ddd = subl_dd_above_ind + ddfract*mass*subld
    subl_dd = subl_dd_above_in + ddfract*subl*mass
    subld = subld - ddfract*subld
    subl = subl - ddfract*subl
    qvd = qvd + evapd + subld
    qv = qv + evap + subl
    ted = ted - cons_alhl*evapd/cons_cp - cons_alhs*subld/cons_cp
    te = te - evap*cons_alhl/cons_cp - subl*cons_alhs/cons_cp
    pfld = mass*qpld
    pfl = qpl*mass
    pfid = mass*qpid
    pfi = qpi*mass
  END IF
  evapd = (qddf3d*evap_dd+qddf3*evap_ddd)/mass
  evap = qddf3*evap_dd/mass
  subld = (qddf3d*subl_dd+qddf3*subl_ddd)/mass
  subl = qddf3*subl_dd/mass
  qvd = qvd + evapd + subld
  qv = qv + evap + subl
  ted = ted - cons_alhl*evapd/cons_cp - cons_alhs*subld/cons_cp
  te = te - evap*cons_alhl/cons_cp - subl*cons_alhs/cons_cp
  qpi = 0.
  qpl = 0.
  pfl_above_outd = pfld
  pfl_above_out = pfl
  pfi_above_outd = pfid
  pfi_above_out = pfi
  evap_dd_above_outd = evap_ddd
  evap_dd_above_out = evap_dd
  subl_dd_above_outd = subl_ddd
  subl_dd_above_out = subl_dd
END SUBROUTINE precipandevap_d

!  Differentiation of marshpalm in forward (tangent) mode:
!   variations   of useful results: diam3 w ve
!   with respect to varying inputs: rain
SUBROUTINE marshpalm_d(rain, raind, pr, diam3, diam3d, ntotal, w, wd, ve&
& , ved)
  IMPLICIT NONE
!Inputs
! in kg m^-2 s^-1, mbar
  REAL*8, INTENT(IN) :: rain, pr
  REAL*8, INTENT(IN) :: raind
!Outputs
  REAL*8, INTENT(OUT) :: diam3, ntotal, w, ve
  REAL*8, INTENT(OUT) :: diam3d, wd, ved
!Locals
  INTEGER :: iqd
!cm^-3
  REAL*8, PARAMETER :: n0=0.08
  REAL*8 :: rain_day, slopr, diam1
  REAL*8 :: rain_dayd
  REAL*8 :: rx(8), d3x(8)
  REAL*8 :: rxd(8), d3xd(8)
  INTRINSIC sqrt
  INTRINSIC max
  REAL*8 :: result1
  rain_day = 0.0
  slopr = 0.0
  diam1 = 0.0
!Marshall-Palmer sizes at different rain-rates: avg(D^3)
!RX = (/ 0.   , 5.   , 20.  , 80.  , 320. , 1280., 5120., 20480. /)  ! rain per in mm/day
  rxd(1) = 0.0_8
  rx(1) = 0.
  rxd(2) = 0.0_8
  rx(2) = 5.
  rxd(3) = 0.0_8
  rx(3) = 20.
  rxd(4) = 0.0_8
  rx(4) = 80.
  rxd(5) = 0.0_8
  rx(5) = 320.
  rxd(6) = 0.0_8
  rx(6) = 1280.
  rxd(7) = 0.0_8
  rx(7) = 5120.
  rxd(8) = 0.0_8
  rx(8) = 20480.
!D3X= (/ 0.019, 0.032, 0.043, 0.057, 0.076, 0.102, 0.137, 0.183  /)
  d3xd(1) = 0.0_8
  d3x(1) = 0.019
  d3xd(2) = 0.0_8
  d3x(2) = 0.032
  d3xd(3) = 0.0_8
  d3x(3) = 0.043
  d3xd(4) = 0.0_8
  d3x(4) = 0.057
  d3xd(5) = 0.0_8
  d3x(5) = 0.076
  d3xd(6) = 0.0_8
  d3x(6) = 0.102
  d3xd(7) = 0.0_8
  d3x(7) = 0.137
  d3xd(8) = 0.0_8
  d3x(8) = 0.183
  rain_dayd = 3600.*24.*raind
  rain_day = rain*3600.*24.
  IF (rain_day .LE. 0.00) THEN
    diam1 = 0.00
    diam3 = 0.00
    ntotal = 0.00
    w = 0.00
    diam3d = 0.0_8
  ELSE
    diam3d = 0.0_8
  END IF
  DO iqd=1,7
    IF (rain_day .LE. rx(iqd+1) .AND. rain_day .GT. rx(iqd)) THEN
      slopr = (d3x(iqd+1)-d3x(iqd))/(rx(iqd+1)-rx(iqd))
      diam3d = slopr*rain_dayd
      diam3 = d3x(iqd) + (rain_day-rx(iqd))*slopr
    END IF
  END DO
  IF (rain_day .GE. rx(8)) THEN
    diam3 = d3x(8)
    diam3d = 0.0_8
  END IF
  ntotal = 0.019*diam3
  diam3d = 0.664*diam3d
  diam3 = 0.664*diam3
  result1 = sqrt(1000./pr)
  wd = result1*2483.8*diam3d
  w = (2483.8*diam3+80.)*result1
  IF (0.99*w/100. .LT. 1.000) THEN
    ve = 1.000
    ved = 0.0_8
  ELSE
    ved = 0.99*wd/100.
    ve = 0.99*w/100.
  END IF
  diam1 = 3.0*diam3
  diam1 = diam1/100.
  diam3d = diam3d/100.
  diam3 = diam3/100.
  wd = wd/100.
  w = w/100.
  ntotal = ntotal*1.0e6
END SUBROUTINE marshpalm_d

!  Differentiation of dqsat_bac in forward (tangent) mode:
!   variations   of useful results: dqsi qssi
!   with respect to varying inputs: temp dqsi qssi
SUBROUTINE dqsat_bac_d(dqsi, dqsid, qssi, qssid, temp, tempd, plo, im, &
& jm, lm, estblx, cons_h2omw, cons_airmw)
  IMPLICIT NONE
!Inputs
  INTEGER :: im, jm, lm
  REAL*8, DIMENSION(im, jm, lm) :: temp, plo
  REAL*8, DIMENSION(im, jm, lm) :: tempd
  REAL*8 :: estblx(:)
  REAL*8 :: cons_h2omw, cons_airmw
!Outputs
  REAL*8, DIMENSION(im, jm, lm) :: dqsi, qssi
  REAL*8, DIMENSION(im, jm, lm) :: dqsid, qssid
!Locals
  REAL*8, PARAMETER :: max_mixing_ratio=1.0
  REAL*8 :: esfac
  INTEGER :: i, j, k
  REAL*8 :: tl, tt, ti, dqsat, qsat, dqq, qq, pl, pp, dd
  REAL*8 :: tld, ttd, tid, dqsatd, qsatd, qqd, ddd
  INTEGER :: it
  INTEGER, PARAMETER :: degsubs=100
  REAL*8, PARAMETER :: tmintbl=150.0, tmaxtbl=333.0
  INTEGER, PARAMETER :: tablesize=nint(tmaxtbl-tmintbl)*degsubs+1
  INTRINSIC nint
  INTRINSIC int
  esfac = cons_h2omw/cons_airmw
  DO k=1,lm
    DO j=1,jm
      DO i=1,im
        tld = tempd(i, j, k)
        tl = temp(i, j, k)
        pl = plo(i, j, k)
        pp = pl*100.0
        IF (tl .LE. tmintbl) THEN
          ti = tmintbl
          tid = 0.0_8
        ELSE IF (tl .GE. tmaxtbl - .001) THEN
          tid = 0.0_8
          ti = tmaxtbl - .001
        ELSE
          tid = tld
          ti = tl
        END IF
        ttd = degsubs*tid
        tt = (ti-tmintbl)*degsubs + 1
        it = int(tt)
        dqq = estblx(it+1) - estblx(it)
        qqd = dqq*ttd
        qq = (tt-it)*dqq + estblx(it)
        IF (pp .LE. qq) THEN
          qsat = max_mixing_ratio
          dqsat = 0.0
          qsatd = 0.0_8
          dqsatd = 0.0_8
        ELSE
          ddd = -((-((1.0-esfac)*qqd))/(pp-(1.0-esfac)*qq)**2)
          dd = 1.0/(pp-(1.0-esfac)*qq)
          qsatd = esfac*(qqd*dd+qq*ddd)
          qsat = esfac*qq*dd
          dqsatd = esfac*degsubs*dqq*pp*(ddd*dd+dd*ddd)
          dqsat = esfac*degsubs*dqq*pp*(dd*dd)
        END IF
        dqsid(i, j, k) = dqsatd
        dqsi(i, j, k) = dqsat
        qssid(i, j, k) = qsatd
        qssi(i, j, k) = qsat
      END DO
    END DO
  END DO
END SUBROUTINE dqsat_bac_d

!  Differentiation of dqsats_bac in forward (tangent) mode:
!   variations   of useful results: dqsi qssi
!   with respect to varying inputs: temp
SUBROUTINE dqsats_bac_d(dqsi, dqsid, qssi, qssid, temp, tempd, plo, &
& estblx, cons_h2omw, cons_airmw)
  IMPLICIT NONE
!Inputs
  REAL*8 :: temp, plo
  REAL*8 :: tempd
  REAL*8 :: estblx(:)
  REAL*8 :: cons_h2omw, cons_airmw
!Outputs
  REAL*8 :: dqsi, qssi
  REAL*8 :: dqsid, qssid
!Locals
  REAL*8, PARAMETER :: max_mixing_ratio=1.0
  REAL*8 :: esfac
  REAL*8 :: tl, tt, ti, dqsat, qsat, dqq, qq, pl, pp, dd
  REAL*8 :: tld, ttd, tid, dqsatd, qsatd, qqd, ddd
  INTEGER :: it
  INTEGER, PARAMETER :: degsubs=100
  REAL*8, PARAMETER :: tmintbl=150.0, tmaxtbl=333.0
  INTEGER, PARAMETER :: tablesize=nint(tmaxtbl-tmintbl)*degsubs+1
  INTRINSIC nint
  INTRINSIC int
  esfac = cons_h2omw/cons_airmw
  tld = tempd
  tl = temp
  pl = plo
  pp = pl*100.0
  IF (tl .LE. tmintbl) THEN
    ti = tmintbl
    tid = 0.0_8
  ELSE IF (tl .GE. tmaxtbl - .001) THEN
    ti = tmaxtbl - .001
    tid = 0.0_8
  ELSE
    tid = tld
    ti = tl
  END IF
  ttd = degsubs*tid
  tt = (ti-tmintbl)*degsubs + 1
  it = int(tt)
  dqq = estblx(it+1) - estblx(it)
  qqd = dqq*ttd
  qq = (tt-it)*dqq + estblx(it)
  IF (pp .LE. qq) THEN
    qsat = max_mixing_ratio
    dqsat = 0.0
    qsatd = 0.0_8
    dqsatd = 0.0_8
  ELSE
    ddd = -((-((1.0-esfac)*qqd))/(pp-(1.0-esfac)*qq)**2)
    dd = 1.0/(pp-(1.0-esfac)*qq)
    qsatd = esfac*(qqd*dd+qq*ddd)
    qsat = esfac*qq*dd
    dqsatd = esfac*degsubs*dqq*pp*(ddd*dd+dd*ddd)
    dqsat = esfac*degsubs*dqq*pp*(dd*dd)
  END IF
  dqsid = dqsatd
  dqsi = dqsat
  qssid = qsatd
  qssi = qsat
END SUBROUTINE dqsats_bac_d

END MODULE cloud_tl