!------------------------------------------------------------------------
!
!  module cape_help
!
!  function hmoist:
!  function enthalpy
!  function rsat
!  subroutine up
!  subroutine t_from_hm
!  subroutine partition
!  function dksdt
!  function drssdt
!
!------------------------------------------------------------------------
module help
!------------------------------------------------------------------------
!
!  use statements
!
!------------------------------------------------------------------------
implicit none
!------------------------------------------------------------------------
!
!  Public Parameters (Thermodynamic Constants)
!
!------------------------------------------------------------------------
double precision, parameter :: p_surf = 1000.     ! assumed surface pressure
double precision, parameter :: rvv = 461.50464    ! should be same as cam
double precision, parameter :: rdd = 287.04231    ! should be same as cam
double precision, parameter :: eps = rdd/rvv
double precision, parameter :: epsi = 1./eps
double precision, parameter :: cpd = 1004.64      ! same as cam Sept 2010
double precision, parameter :: cvd = 717.
double precision, parameter :: cl = 4188.0        ! same as cam4
double precision, parameter :: cpv = 1810.0       ! same as cam Sept 2010
double precision, parameter :: ci = 2117.27       ! same as cam4
double precision, parameter :: cpvmcl = cpv - cl
double precision, parameter :: clmci = cl - ci
double precision, parameter :: lv0 = 2.501E+06    ! SHR_CONST_LATVAP in cam4
double precision, parameter :: lf0 = 3.337E+05    ! SHR_CONST_LATICE in cam4
double precision, parameter :: tkelvin = 273.15
double precision, parameter :: g = 9.80616         ! same as cam Sept 2010
double precision, parameter :: ginv = 1.0/g
double precision, parameter :: p0 = 100000.   ! cam4 has: SHR_CONST_PSTD    = 101325.0_R8  standard pressure
double precision, parameter :: rhow = 1000.   ! density of water
double precision, parameter :: bad = -999.
!------------------------------------------------------------------------
!
!   Allowable errors
!   - these are also used in conv_tend
!   - In double precision, errors can be set at least 100 times smaller
!   - these should all be relative errors
!
!  rherr_max: only used in subroutine t_from_hm. Max allowed deviation of rv from
!     rvs when look for saturated solution. Error here does not give rise to
!     errors in total water
!
!  Reduced hmerr_max to 1.0E-07 due to problems in subroutine start. It is likely
!   that these problems are due to residual discontinuity in esat at 0 C, if this
!   was fixed could probably lower the allowed error again.
!
!  Sep 2010: increased rherr_max to 0.01. Since was generaing an error. This
!     is not a real "error", so shouldn't really matter. Just strange that parcels
!     could be so supersaturated.
!
!------------------------------------------------------------------------
double precision, parameter :: hmerr_max = 1.0E-06   ! max hm/km error for double precision
double precision, parameter :: rterr_max = 1.0E-06   ! May 2010: inc for starting spectrum (FIX)
double precision, parameter :: masserr_max = 1.0E-09 ! max rh error for double precision
double precision, parameter :: rherr_max = 1.0E-01   ! max rh error; only used in subroutine t_from_km
double precision, parameter :: dperr_max = 1.0E-09   ! max dp error in tendency calc
!------------------------------------------------------------------------
!
!   Public interfaces
!
!------------------------------------------------------------------------
  public t_from_km ! (subroutine) 
  public partition ! (subroutine)
  public esat      ! (function)
  public rsat      ! (function)
  public drssdt    ! (function)
  public enthalpy  ! (function)
  public hmoist    ! (function)
  public latent    ! (function)
  public fusion    ! (function)
  public dksdt     ! (function)
!------------------------------------------------------------------------
!
!   contains
!
!------------------------------------------------------------------------
contains
!------------------------------------------------------------------------
!**********************************************************************73
!------------------------------------------------------------------------
!
!  function hmoist:
!
!  Inputs: t (K), rt, rv, ri, z (m)
!
!  hm = (cpd + rt*cl)*T + lv*rv - lf*ri + (1. + rt)*g*z
!
!  [hm] = [J/kg dry air]
!
!  Note that cv does not enter explicitly but only implicitly through lv.
!
!  The input z should be the mass weighted z of a layer. Simple mean height
!  could introduce errors.
!
!------------------------------------------------------------------------
double precision function hmoist(tt,rt,rv,ri,zz)
implicit none
double precision, intent(in) :: tt,rt,rv,ri,zz
  hmoist = enthalpy(tt,rt,rv,ri) + (1. + rt)*g*zz
end function hmoist
!------------------------------------------------------------------------
!**********************************************************************73
!------------------------------------------------------------------------
!
!   function sigmoidal
!
!------------------------------------------------------------------------
double precision function sigmoidal(x_in,x_half,x_scale,y_min,y_max)
!------------------------------------------------------------------------
!
!   Use Statements
!
!------------------------------------------------------------------------
implicit none
!------------------------------------------------------------------------
!
!   In/Out variables
!
!------------------------------------------------------------------------
double precision, intent(in) :: x_in
double precision, intent(in) :: x_half
double precision, intent(in) :: x_scale
double precision, intent(in) :: y_min
double precision, intent(in) :: y_max
!------------------------------------------------------------------------
!
!
!------------------------------------------------------------------------
double precision :: tt
!------------------------------------------------------------------------
!
!
!------------------------------------------------------------------------
tt = (x_in - x_half)/x_scale
sigmoidal = y_min + (y_max/(1. + exp(-tt)))
!------------------------------------------------------------------------
!
!    END     (function sigmoidal)
!
!------------------------------------------------------------------------
end function sigmoidal
!------------------------------------------------------------------------
!**********************************************************************73
!------------------------------------------------------------------------
!
!  Function enthalpy (uses definition 4.9.2 of Emanuel, includes ice)
!
!  k = (cpd + rt*cl)*T + lv*rv - lf*ri
!
!  Inputs: t (K), rt, rv, ri
!
!------------------------------------------------------------------------
double precision function enthalpy(tt,rt,rv,ri)
implicit none
double precision, intent(in) :: tt,rt,rv,ri
double precision :: lv,lf
  lv = latent(tt)
  lf = fusion(tt)
  enthalpy = (cpd + rt*cl)*tt + lv*rv - lf*ri
end function enthalpy
!------------------------------------------------------------------------
!**********************************************************************73
!------------------------------------------------------------------------
!
!  function rsat
!
!  Inputs:
!    t: temperature (K)
!    p: pressure (Pa)
!
!  Outputs:
!    rsat
!
!  Interesting: at low pressures in the stratosphere, es > p,
!  so that rsat becomes negative. In this case define rsat as eps*es/p
!
!  Don't think rsat should go negative at ordinary range of t/p
!
!------------------------------------------------------------------------
double precision function rsat(tt,pp)
implicit none
double precision, intent(in) :: tt,pp
double precision :: es
  es = esat(tt)
  if (pp > 500.) then
    rsat = eps*es/(pp - es)
  else
    rsat = eps*es/pp
  endif
  if (rsat < 0.) then
     print*,'======== rsat negative in function rsat ==========================='
     print*,'sometimes a result of extremely high unrealistic temperatures of around 370 K'
     print*,'tt pp = ',tt,pp
     print*,'rsat = ',rsat
     print*,'eps es pp =',eps,es,pp
     stop 'rsat gone negative'
  endif
end function rsat
!------------------------------------------------------------------------
!**********************************************************************73
!------------------------------------------------------------------------
!
!  kmm = (cpd + rtm*cl)*tm + lv*rvm - lf*rim
!  kmm = (cpd + rvm*cl)*tm + lv*rvm   (rim = rlm = 0.)
!  kmm - lv*rvm = (cpd + rvm*cl)*tm   (rim = rlm = 0.)
!  kmm - (lv0 + cpvmcl*tc)*rvm = (cpd + rvm*cl)*tm   
!  kmm - (lv0 + cpvmcl*(tm - tkelvin))*rvm = (cpd + rvm*cl)*tm   
!  kmm - (lv0 + cpvmcl*tm - cpvmcl*tkelvin)*rvm = (cpd + rvm*cl)*tm   
!  kmm - lv0*rvm - cpvmcl*tm*rvm + cpvmcl*tkelvin*rvm = (cpd + rvm*cl)*tm   
!  kmm - lv0*rvm + cpvmcl*tkelvin*rvm = (cpd + rvm*cl)*tm + cpvmcl*tm*rvm
!  kmm + rvm*(cpvmcl*tkelvin - lv0) = tm*(cpd + rvm*cl + cpvmcl*rvm)
!  kmm + rvm*(cpvmcl*tkelvin - lv0) = tm*(cpd + rvm*cl + cpvmcl*rvm)
!
!------------------------------------------------------------------------
subroutine t_from_km_allow_sat(kmm,rvm,pb,tm)
!------------------------------------------------------------------------
!
!  In/out variables
!
!------------------------------------------------------------------------
implicit none
double precision, intent(in) :: kmm,rvm,pb
double precision, intent(out) :: tm
!------------------------------------------------------------------------
!
!  Local variables
!
!------------------------------------------------------------------------
double precision :: AA,BB
!------------------------------------------------------------------------
!
!   Determine solution
!
!------------------------------------------------------------------------
AA = kmm - rvm*(lv0 - cpvmcl*tkelvin)
BB = cpd + rvm*cl + cpvmcl*rvm
tm = AA/BB
!------------------------------------------------------------------------
!
!------------------------------------------------------------------------
end subroutine t_from_km_allow_sat
!------------------------------------------------------------------------
!**********************************************************************73
!------------------------------------------------------------------------
!
!  subroutine t_from_km
!
!  Inputs:
!  (i) kmm: total enthalpy of parcel
!  (ii) rtm: total water of parcel
!  (iii) pb: background pressure
!
!  Outputs:
!  (i) rvm
!  (ii) rlm
!  (iii) rim
!  (iv) tm
!
!
!  Four Constraints:
!  (i) total water conservation
!  (ii) moist enthalpy conservation
!       km = (cpd + cl*rt)*T + rv*Lv - ri*Lf (Emanuel definition, page 128)
!  (iii) C-C relation
!  (iv) assumed ice/(ice + liquid water) partitioning.
!       The ri/(rl + ri) ratio is an externally imposed, temperature
!       dependent constraint.
!
!  The saturation vapor pressure does not depend on the presence of
!  ice, i.e. es = esi for T < 0 C, and es = esw for T > 0 C. This
!  is inconsistent with how ice is treated, i.e. allowing for the existence
!  of supercooled water. However, allowing es to depend on the presence
!  of ice would introduce a discontinuity in es which would likely make
!  the numerical solution much more difficult.
!
!------------------------------------------------------------------------
subroutine t_from_km(kmm,rtm,pb,rvm,rlm,rim,tm)
!------------------------------------------------------------------------
!
!  use statements
!
!------------------------------------------------------------------------
!------------------------------------------------------------------------
!
!  In/out variables
!
!------------------------------------------------------------------------
implicit none
double precision, intent(in) :: kmm,rtm,pb
double precision, intent(out) :: tm
double precision, intent(inout) :: rvm,rlm,rim
!------------------------------------------------------------------------
!
!  Local variables
!
!------------------------------------------------------------------------
double precision :: tg,kg,terror,kmerror,aerror,dkgdt,dk,khave,rsg,rvg,rlg,rig
double precision :: rherror,rhmax,dt,rsm,rterror
double precision :: rtt,rterr,small,AA,BB,tstart,rl_in,rv_in
integer :: niter_max,iter
!------------------------------------------------------------------------
!
!  Store input variables for diagnostics
!
!------------------------------------------------------------------------
rv_in = rvm
rl_in = rlm
!print*,'kmm = ',kmm
!print*,'rvm = ',rvm
!print*,'rlm = ',rlm
!print*,'rim = ',rim
!print*,'rtm = ',rtm
!print*,'tm = ',tm
!print*,'pb = ',pb
!------------------------------------------------------------------------
!
!  Check for negative kmm or rtm on entry:
!
!------------------------------------------------------------------------
  if (kmm < 0.) then
    stop 'negative kmm on entry to subroutine t_from_km'
  endif
  if (rtm < 0.) then
    print*,'WARNING: negative rtm on entry to subroutine t_from_km'
    print*,'rtm = ',rtm
    print*,'input rvm = ',rvm
    print*,'input rlm = ',rlm
!    print*,'setting rtm equal to zero'
!    rtm = 0.
    stop 'negative rtm on entry to subroutine t_from_km'
  endif
  if ((rvm+rim+rlm)<= 0.) then
    print*,'rvm = ',rvm
    print*,'rim = ',rim
    print*,'rlm = ',rlm
    print*,'WARNING: negative (rvm+rim+rlm) on entry to subroutine t_from_km'
    stop 'negative (rvm+rim+rlm) on entry to subroutine t_from_km'
  endif
!------------------------------------------------------------------------
!
!  First: check to make sure input rvm/rlm/rim consistent with input rtm
!
!  Calculate a temperature tm from enthalpy kmm assumming fixed condensate
!  (i.e. rtm not equal to rvm)
!
!  kmm = (cpd + rtm*cl)*tm + lv*rvm - lf*rim
!  (cpd + rtm*cl)*tm = kmm - lv*rvm + lf*rim
!  kmm - [lv0 - cpvmcl*tkelvin]*rvm + [lf0 - clmci*tkelvin]*rim
!  tm*(cpd + rtm*cl + cpvmcl*rv - clmci*rim) =
!  kmm - rvm*(lv0 - cpvmcl*tkelvin) + rim*(lf0 - clmci*tkelvin)
!  Gives:
!  tm = AA/BB, where:
!  AA = kmm - rvm*(lv0 - cpvmcl*tkelvin) + rim*(lf0 - clmci*tkelvin)
!  BB = (cpd + rtm*cl + cpvmcl*rv - clmci*rim)
!
!------------------------------------------------------------------------
!------------------------------------------------------------------------
!
!  rtm must be defined on entry
!  rvm/rlm/rim in general undefined: DO NOT USE
!
!  Look for tm where all condensate evaporated
!  Assume rtm = rvm, rim = 0.
!
!------------------------------------------------------------------------
  rvm = rtm
  AA = kmm - rvm*(lv0 - cpvmcl*tkelvin)
  BB = cpd + rvm*cl + cpvmcl*rvm
  tm = AA/BB
!------------------------------------------------------------------------
!
!  Have a guess for tm:
!  - from tm -> rsm
!  - check rvm < rs(tm).
!  - If yes, have solution
!
!------------------------------------------------------------------------
  rsm = rsat(tm,pb)
  if (rvm <= rsm) then
!    print*,'unsaturated solution in t_from_km works'
!------------------------------------------------------------------------
!
!   Check error:
!
!------------------------------------------------------------------------
    khave = enthalpy(tm,rtm,rvm,rim)
    kmerror = abs(khave-kmm)/kmm
    if (kmerror > hmerr_max) then
      print*,'----------- large kmerror in t_from_km ----------------'
      print*,'relative error in enthalpy = ',kmerror
      print*,'kmm = ',kmm
      print*,'khave = ',khave
      print*,'This is for the unsaturated solution'
      print*,'in unsaturated solution'
      stop 'relative error in enthalpy in t_from_km too big'
    endif
  else
!------------------------------------------------------------------------
!
!   Saturated Solution (rtm > rs(tm)):
!   ---------------------------------
!
!   Keep rtm and kmm fixed
!   Let tg = tm (unsaturated guess)
!   call partition: assume tg is true and use pb + rtm to find rvg,rlg,rig,rsg
!   find kg
!   loop over iter:
!     Let tm = tp + dT
!     dk = kmm - kg
!     define dkgdt (assumes parcel saturated)
!     dt = dk*(dT/dk)
!     tg = tg + dt (corrected tg)
!     find new kg
!     get estimate of t error, relative k error
!    end loop
!
!  Should exit loop with reasonable estimate of tg
!
!  Sept 2009: was having problems with rh > 1.002 on exit from t_from_km, and
!             tried to fix by increasing niter from 5 to 7, but didn't
!             seem to help.
!
!  Can't really have a while loop here since no real error test.
!
!------------------------------------------------------------------------
    niter_max = 6
    iter = 0
    tg = tm
!    print*,'calling partition from t_from_km'
    call partition(tg,pb,rtm,rvg,rlg,rig,rsg)
    kg = enthalpy(tg,rtm,rvg,rig)
    kmerror = 1.
    do iter = 1,niter_max
      dk = kmm - kg
      dkgdt = dksdt(tg,pb,rtm,rsg,rig)
      dt = dk/dkgdt
      tg = tg + dt
      call partition(tg,pb,rtm,rvg,rlg,rig,rsg)
      kg = enthalpy(tg,rtm,rvg,rig)
!------------------------------------------------------
!  rterror - should be zero for every iteration
!  the rvg rlg rig that come out of partition should sum to rtm
!------------------------------------------------------
      rterror = abs((rvg + rlg + rig) - rtm)/rtm
      if (rterror > 1.0E-06) then
         print*,'---------- rterror in t_from_km should be zero for every iteration -----'
         print*,'rterror in t_from_km = ',rterror
         print*,'iter = ',iter
         print*,'niter_max = ',niter_max
         print*,'entry water variables : rtm rvm rl rim'
         print*,'rtm = ',rtm
         print*,'rvm = ',rvm
         print*,'rlm = ',rlm
         print*,'rim = ',rim
         print*,'rvg = ',rvg
         print*,'rlg = ',rlg
         print*,'rig = ',rig
         print*,'tg = ',tg
         print*,'dk = ',dk
         print*,'kmm = ',kmm
         print*,'kg = ',kg
         print*,'dkgdt = ',dkgdt
         print*,'------------------------------------------------------------'
         stop 'rterror in subroutine t_from_km during iteration'
      endif
    end do
!------------------------------------------------------------------------
!
!  Want a solution which preserves total water and moist enthalpy.
!
!  If believe tg, will be some error in moist enthalpy
!  Instead, believe rtm partitioning: rvg,rlg,rig,rsg.
!  This guarantees water conservation.
!  Find tm by enforcing km conservation.
!  Error will be that rs(tm) will not exactly equal rvg.
!  Just have to live with the possibility of small super/under saturations
!   for air parcels that have condensate.
!
!  km = (cpd + cl*rtm)*tm + rvm*lv(tm) - rim*lf(tm)
!  km = (cpd + cl*rtm)*tm + rvm*(lv0 + cpvmcl*tc) - rim*(lf0 + clmci*tc)
!  km = (cpd + cl*rtm)*tm + rvm*[lv0 + cpvmcl*(tm - tkelvin)]
!       - rim*[lf0 + clmci*(tm - tkelvin)]
!  km = (cpd + cl*rtm + rvm*cpvmcl - rim*clmci)*tm + rvm*[lv0 - cpvmcl*tkelvin]
!       - rim*[lf0 - clmci*tkelvin]
!  km - rvm*[lv0 - cpvmcl*tkelvin] + rim*[lf0 - clmci*tkelvin] =
!       (cpd + cl*rtm + rvm*cpvmcl - rim*clmci)*tm
!  tm = [km + rvm*(cpvmcl*tkelvin - lv0) + rim*(lf0 - clmci*tkelvin)]/
!       (cpd + cl*rtm + rvm*cpvmcl - rim*clmci)
!
!  This expression for tm should exactly conserve km. Check that rterror = 0.
!  The problem is that rsm from tm will not exactly equal rvg, even though
!  the parcel is assumed saturated.
!  Check rherror
!
!------------------------------------------------------------------------
  rvm = rvg
  rim = rig
  rlm = rtm - rvm - rim
! print*,'rim = ',rim
! print*,'rtm = ',rtm
! print*,'rvm = ',rvm
! print*,'rim = ',rim
!------------------------------------------------------------------------
!
!  Check if relative error in rtp is OK (should be zero)
!
!------------------------------------------------------------------------
  rterror = abs((rvm + rlm + rim) - rtm)/rtm
!  print*,'final rterror in t_from_km = ',rterror
  if (rterror > 1.0E-06) then
    print*,'------------- rterror in t_from_km should be zero ----------------'
    print*,'This is for the saturated solution'
    print*,'rterror = ',rterror
    print*,'rherror = ',rherror
    print*,'input rtm = ',rtm
    print*,'input rv_in = ',rv_in
    print*,'input rl_in = ',rl_in
    print*,'solution rvm = ',rvm
    print*,'solution rlm = ',rlm
    print*,'solution rim = ',rim
    print*,'should be same as rvm: rvg = ',rvg
    print*,'should equal rlm: rtm - rvm - rim = ',rtm - rvm - rim
    print*,'rsm = ',rsm
    print*,'tm = ',tm
    print*,'khave = ',khave
    print*,'-----------------------------------------------------------------'
    stop 'relative rt error in t_from_km non-zero'
  endif
!------------------------------------------------------------------------
!
!  Use saturated rvm to find tm
!  tm = AA/BB, where:
!  AA = kmm - rvm*(lv0 - cpvmcl*tkelvin) + rim*(lf0 - clmci*tkelvin)
!  BB = (cpd + rtm*cl + cpvmcl*rv - clmci*rim)
!
!------------------------------------------------------------------------
  AA = kmm + rvm*(cpvmcl*tkelvin - lv0) + rim*(lf0 - clmci*tkelvin)
  BB = cpd + cl*rtm + rvm*cpvmcl - rim*clmci
!  print*,'AA = ',AA
!  print*,'BB = ',BB
  tm = AA/BB
!  print*,'tm = ',tm
  khave = enthalpy(tm,rtm,rvm,rim)
!  print*,'khave = ',khave
  kmerror = abs(khave-kmm)/khave
!  print*,'kmerror = ',kmerror
  rsm = rsat(tm,pb)
  rherror = abs(rsm - rvm)/rsm
!------------------------------------------------------------------------
!
!   Check if relative error in moist enthalpy OK (should be zero)
!   Designed to exactly conserve enthalpy, so error should be zero.
!   - probably redundant
!
!------------------------------------------------------------------------
!  print*,'final kmerror in t_from_km = ',kmerror
  if (kmerror > hmerr_max) then
    print*,'------- enthalpy error in subroutine t_from_km exceeds tolerance -----'
    print*,'This routine should exactly conserve total water and enthalpy by construction'
    print*,'If conservation is violated might be a segmentation fault'
    print*,'The only error should be that the desired input rh does not equal actual rh'
    print*,'this difference is rherror = ',rherror
    print*,'relative error in enthalpy = ',kmerror
    print*,'This is for the saturated solution'
    print*,'max relative error = ',hmerr_max
    print*,'kmm = ',kmm
    print*,'khave = ',khave
    print*,'rterror = ',rterror
    print*,'input rtm = ',rtm
    print*,'rvm = ',rvm
    print*,'rsm = ',rsm
    print*,'rlm = ',rlm
    print*,'rim = ',rim
    print*,'tm = ',tm
    print*,'AA = ',AA
    print*,'BB = ',BB
    print*,'cpvmcl = ',cpvmcl
    print*,'cl = ',cl
    print*,'cpvmcl times tkelvin - lv0 = ',cpvmcl*tkelvin - lv0
    print*,'---------------------------------------------------------------'
    stop 'relative enthalpy error in subroutine t_from_km too big'
  endif
!------------------------------------------------------------------------
!
!   Check if relative error in rh is small (in general will be non-zero)
!   - Forced to allow some small supersaturation in forcing rt/km conservation.
!   - Program should just convert any supersaturation to rl, so should not cause
!     a problem.
!   - Jan 25/2011: getting rh errors at very cold Temperatures (200 K).
!
!------------------------------------------------------------------------
  if (rherror > rherr_max) then
     print*,'------- WARNING: RH error exceeds threshold in subroutine t_from_km --------'
     print*,'This is for the saturated solution'
     print*,'May not be that serious: can still have small km and rt errors'
     print*,'May mean that original assumption of saturated solution not valid'
     print*,'pressure in hPa = ',pb*0.01
     print*,'tm = ',tm
     print*,'rsm = ',rsm
     print*,'rvm = ',rvm
     print*,'rlm = ',rlm
     print*,'input rtm = ',rtm
     print*,'input rv_in = ',rv_in
     print*,'input rl_in = ',rl_in
     print*,'rherror = abs(rsm - rvm)/rsm'
     print*,'input moist enthalpy kmm = ',kmm
     print*,'final moist enthalpy khave = ',khave
     print*,'relative error in km kmerror = ',kmerror
     print*,'relative error in RH = ',rherror
     print*,'rherr_max = ',rherr_max
     print*,'final relative total water = ',rterror
     print*,'Try to figure out why unsaturated solution was not chosen'
     print*,'----------------------------------------------------------'
!     stop 'relative error in RH in t_from_km too big'
  endif
!------------------------------------------------------------------------
!
!  endif for sat/unsat solution
!
!------------------------------------------------------------------------
  endif
!------------------------------------------------------------------------
!
!   End
!
!------------------------------------------------------------------------
end subroutine t_from_km
!------------------------------------------------------------------------
!**********************************************************************73
!------------------------------------------------------------------------
!
!   Subroutine partition
!
!   Inputs: t (K), p (Pa), rt
!
!   Outputs: rv, rl, ri, rs
!
!  Treats that t/p/rt as fixed.
!  Partitions total water under the assumption that all condensate is
!    evaporated if possible, i.e parcel is saturated if rt large enough.
!  Since the re-partitioning of water that is done will violate km conservation,
!    partition should only be called if km conservation is being enforced in
!    some other way.
!  Should conserve total water.
!
!  Should only be called from t_from_km when it is known that system is saturated
!
!   Nov/2009: revised to not allow ice
!
!   Checks that system is saturated, i.e. rt > rs
!
!------------------------------------------------------------------------
subroutine partition(tt,pp,rt,rv,rl,ri,rs)
implicit none
double precision, intent(in) :: tt,pp,rt
double precision, intent(out) :: rv,rl,ri,rs
double precision :: small,rterror,rtabs
!------------------------------------------------------------------------
!
!
!------------------------------------------------------------------------
  rs = rsat(tt,pp)
  if (rt <= rs) then
    ri = 0.
    rv = rt
    rl = 0.
  else
    ri = 0.
    rv = rs
    rl = rt - rs - ri
  endif
!------------------------------------------------------------------------
!
!  In the formula  rl = rt - ri - rs, rounding errors can make rl < 0.
!  If rl is negative by a very small amount, set to zero.
!
!------------------------------------------------------------------------
  small = 1.0E-08
  if ((rl < 0.).and.(rl > (-1.*small))) then
     rv = rv + rl
     rl = 0.
     print*,'setting rl to zero in partition'
  endif
  if (ri < 0.) stop 'ri problem in partition'
  if (rl < 0.) then
    print*,'tt pp rt rv rs rl ri = ',tt,pp,rt,rv,rs,rl,ri
    stop 'rl problem in partition'
  endif
!------------------------------------------------------------------------
!
!   End
!
!------------------------------------------------------------------------
end subroutine partition
!------------------------------------------------------------------------
!**********************************************************************73
!------------------------------------------------------------------------
!
!   Function dksdt
!
!   Evaluates partial derivative of km with respect to T, assuming the parcel
!    is saturated, and that rt is fixed.
!
!   See notes for the derivation of:
!
!   drs/dT = des/dT*(eps*p)/(p - es)**2
!   des/dT = Lv*es/(Rv*T*T)  (valid only for water?)
!   dkm/dT = cpd + rtm*cl + cpvmcl*rsm + lv*drsmdt - clmci*rim -lf*drimdt
!
!   dk/dT = cpd + rt*cl + cpvmcl*rs + lv*drsdt - clmci*ri -lf*dridt
!
!   Nov 2009: removed dridt term, so now not valid for non-zero ri
!
!------------------------------------------------------------------------
double precision function dksdt(tt,pp,rt,rs,ri)
implicit none
double precision, intent(in) :: tt,pp,rt,rs,ri
double precision :: lv,lf
  lv = latent(tt)
  lf = fusion(tt)
  dksdt = cpd + rt*cl + cpvmcl*rs + lv*drssdt(tt,pp) - clmci*ri
end function dksdt
!------------------------------------------------------------------------
!**********************************************************************73
!------------------------------------------------------------------------
!
!  Function drssdt
!
!  Inputs:
!    t: temperature (K)
!    p: pressure (Pa)
!
!  Outputs:
!    drssdt
!
!  If lf is not included below 0 C, desdt is in error by about 10%. With
!    lf, the error is about 0.1%.
!
!------------------------------------------------------------------------
double precision function drssdt(tt,pp)
implicit none
double precision, intent(in) :: tt,pp
double precision :: lv, lf, es, tc, desdt
double precision :: rsat
!------------------------------------------------------------------------
!
!
!
!------------------------------------------------------------------------
  lv = latent(tt)
  lf = fusion(tt)
  es = esat(tt)
  tc = tt - tkelvin
  if (tc > 0.) then
    desdt = (lv*es)/(rvv*tt*tt)
  else
    desdt = ((lv+lf)*es)/(rvv*tt*tt)
  endif
  drssdt = (desdt*eps*pp)/((pp - es)**2)
end function drssdt
!------------------------------------------------------------------------
!**********************************************************************73
!------------------------------------------------------------------------
!
!   Function esat
!
!   Input: temperature tt (K)
!   Output: saturation vapor pressure in Pa (not mb), with respect to water
!      for T > 0 C, and with respect to ice for T < 0 C.
!   This definition of esat may be inconsistent with the the possible allowance
!     of the existence of supercooled water in the treatment of ice. However,
!     allowing for the possibility of a discontinuity in esat likely would
!     lead to numerical difficulties.
!
!   esat formula taken from voemel web site, WMO, 2008
!   gives 6.112 hPa at 0 C
!
!  original expression for esat gave:
!   esat = exp(23.33086 - (6111.72784/t) + 0.15215*log(t))
!  But this did not exactly give es = 6.112 hPa at 0 C, so changed it a bit.
!  Discontinuity in es at 0 C was causing problems in subroutine start.
!
!  In fortran log mean natural log, so is the log t in the ice esat expression
!    really ln?? may have made a mistake, and caused by rh gap due to a discontinuity
!    in esat at 0 C.
!
!  need: log(e0) = AA - (BB/t) + CC*log(t)
!        AA = log(e0) + (BB/t) - CC*log(t)
!
!------------------------------------------------------------------------
double precision function esat(tt)
implicit none
double precision, intent(in) :: tt
double precision :: tc,denom
  tc = tt - tkelvin
  denom = 243.12 + tc
  if (tc > 0.) then
    esat = 6.112*exp(17.62*tc/denom)
  else
    esat = exp(23.33167801703645 - (6111.72784/tt) + 0.15215*log(tt))
  endif
  esat = 100.*esat
end function esat
!------------------------------------------------------------------------
!**********************************************************************73
!------------------------------------------------------------------------
!
!   Latent Heat of Evaporation
!
!   Input: t (K)
!
!------------------------------------------------------------------------
double precision function latent(tt)
implicit none
double precision, intent(in) :: tt
double precision :: tc
  tc = tt - tkelvin
  latent = lv0 + cpvmcl*tc
end function latent
!------------------------------------------------------------------------
!**********************************************************************73
!------------------------------------------------------------------------
!
!   Latent Heat of Fusion
!
!   Input: tt (K)
!
!   cl = 4188.0
!   ci = 2117.27
!   clmci = 2070.73
!   Sept 2010: changed to set equal to lf0 for T > tkelvin
!   - increasing function of T
!
!------------------------------------------------------------------------
double precision function fusion(tt)
implicit none
double precision, intent(in) :: tt
double precision :: tc
  tc = tt - tkelvin
  fusion = lf0 + clmci*tc
  if (tc > 0.) fusion = lf0
end function fusion
!------------------------------------------------------------------------
!**********************************************************************73
!------------------------------------------------------------------------
!
!   end module
!
!------------------------------------------------------------------------
end module help
!------------------------------------------------------------------------
!**********************************************************************73
!------------------------------------------------------------------------
