! test program
implicit none

  real :: Ta, RH, Pa
  real :: q			! specific humidity (kg/kg)
  real :: Rho_v			! vapor density (kg/m3)
  real :: Tdew 			! dew point temp (C. deg)
  real :: e			! water vapor press. in hPa

  real :: Tw                    ! wet bulb temperature (C. deg)

 integer :: i

  Ta = 5.0			! C. deg.
  RH = 90.0			! 0-100 %
  Pa = 1000.0			! hPa


!  call wetbulbtemp(Tw, Ta, RH, Pa)
!  write(*,*) 'Tw = ', Tw
  call wetbulbtemp(Tw, Ta, RH, Pa)
  write(*,*) 'Tw = ', Tw

! do i=1, 60
!  Ta = 50-i
!  call specific_hum( Ta, 100.0, Pa, q )
!  write(*,*) Ta, q
! end do
  


!  call satvap2(Ta, e)
!  e = e * (RH/100.0)
!  call calc_Tdew( e, Tdew )
 ! write(*,*) "vapor press(hPa) =", e, "Dew point Temp.(C.deg)", Tdew

!  call specific_hum( Ta, RH, Pa, q )

!  write(*,*) 'Temp =', Ta, 'C.deg.', 'RH = ', RH, '%',&
!  	'Pa = ', Pa, '(hPa)'
!  write(*,*) 'Specific humidity = ', q , ' (kg/kg)'


! call vapor_density( Ta, RH, Rho_v )
!  write(*,*) 'Rho_v=', Rho_v, 'kg/m3'

!  call sphum2rhov( q, Ta, Pa, Rho_v )
!  write(*,*) 'Rho_v=', Rho_v, 'kg/m3'

  stop
  end



subroutine wetbulbtemp(Tw, Ta, RH, Pa)
! 2015 12 16 
implicit none
 real :: Tw, Ta, RH, Pa ! Tw, Ta in C deg. RH in %, Pa hPa

 real :: dt
 real :: q
! real, parameter :: Cpd=1004.0, Cpv=1847.0, T_ZERO = 273.15
! real, parameter ::  Latent_evap=2.5e6, Latent_melt = 0.3335e6
 real, parameter :: lim = 0.001
 real ::  dif1, dif2, T1, T2, T
 

  call specific_hum( Ta, RH, Pa, q )

  dt = 4.0
  T1 = Ta; T2 = Ta - 0.5*dt

do while( abs(dt) >  lim )
! write(*,*) T1, T2
  call diff_wetbulb_func( q, Ta, Pa, T1, dif1 )
  call diff_wetbulb_func( q, Ta, Pa, T2, dif2 )

 if(dif2 == 0.0 ) then
   Tw = T2
   exit
 else if( dif1 == 0.0 ) then
   Tw = T1
   exit
 else
  if( dif1* dif2 > 0 ) then
     if( abs(dif1) >= abs(dif2) ) then
       T1 = 2*T2 - T1
       dt = abs(T1-T2)

     else!   if( abs(dif2) > abs(dif1) ) 
       T2 = 2*T1 - T2
       dt = abs(T1-T2)
     end if
  else  ! different sign 
     if( abs(dif1) >= abs(dif2) ) then
       T1 = (T2+T1)/2.0
       dt = (T2-T1)
       if( abs(dt) < lim ) then
        Tw = T1
       exit
       end if
     else
       T2 = (T2+T1)/2.0
       dt = (T1-T2)
       if( abs(dt) < lim ) then
        Tw = T2
       exit
       end if
     end if
 end if
end if
end do

end

subroutine diff_wetbulb_func( q, Ta, Pa, T, dif )
implicit none
 real :: q, Ta, Pa, dif, T
 real, parameter :: Cpd=1004.0, Cpv=1847.0, T_ZERO = 273.15
 real, parameter ::  Latent_evap=2.5e6, Latent_melt = 0.3335e6
 real :: le1, h1   ! (LE-H)/(\rho U Ch) = L(qsat(Tw)-qa) - Cp(Ta-Tw)
 real :: Cp, qsat, Latent

 Cp = Cpd*(1-q) + Cpv*q
 if( T >= 0.0 ) then
  Latent = Latent_evap
 else
  Latent = Latent_evap + Latent_melt
 end if
  call specific_hum( T, 100.0, Pa, qsat )

 le1 =  Latent * ( qsat - q )
 h1 =   Cp * (Ta- T )
 dif = le1 - h1
!  write(*,*) "DIF", dif,":", le1, h1
end 


! subroutine wetbulbtemp(Tw, Ta, RH, Pa)
! implicit none
!  real :: Tw, Ta, RH, Pa ! Tw, Ta in C deg. RH in %, Pa hPa
! 
!  real :: tmp, dt
!  real :: q, Cp, Latent, qsat
!  real, parameter :: Cpd=1004.0, Cpv=1847.0, T_ZERO = 273.15, C_ZERO=0
!  real, parameter ::  Latent_evap=2.5e6, Latent_melt = 0.3335e6
!  real, parameter :: lim = 0.01
!  real :: dqsatdt
!  
! 
!   call specific_hum( Ta, RH, Pa, q )
!   tmp = Ta; dt = 1
! 
! do while( abs(dt) >  lim )
!   Cp = Cpd*(1-q)+Cpv*q
!   if( tmp > C_ZERO ) then
!    Latent = Latent_evap
!   else
!    ! Latent = Latent_evap
!    Latent = Latent_evap + Latent_melt
!   end if
! 
!   call specific_hum( tmp, 100.0, Pa, qsat )
!   dt = ((Cp/Latent)*(Ta-tmp) + qsat -q)/(Cp/Latent - dqsatdt(tmp, Pa))
!   tmp = tmp + dt
!  write(*,*) Ta,RH,tmp,dt
! end do
! 
! Tw = tmp
! end subroutine 
! 


subroutine vapor_density( Ta, RH, Rho_v )
  implicit none
  real :: Ta, RH, Rho_v
  real, parameter :: Rd=287.0, Rw = 461.0
  real :: p, tk, e, rh1

  tk = Ta + 273.15
  rh1 = 0.01*RH

  call satvap2(Ta, e) ; e = 100.0*e*rh1

  Rho_v = e/(Rw*tk)
end


subroutine sphum2rhov(  q, Ta, Pa, Rho_v )
!  calc.  RH from q (kg/kg), Ta (C deg.), Pa (hPa)
!
!   p - e = rho_a Rd T
!       e = rho_v Rw T
!   p = rho_a Rd T + rho_v Rw T
!   q = rho_v/( rho_a + rho_v)
!     rho_a = q/(1-q) * rho_v
!   rho_v = p/(T  ( (1-q)/q * Rd + Rw ) )


   implicit none
   real ::  q, Ta, Pa, RH
   real ::  Rho_v
   real :: p, tk, q1
   real, parameter :: Rd=287.0, Rw = 461.0

   p = Pa * 100.0;  tk = Ta + 273.15
   q1 = (1.0-q)/q

   Rho_v = p/( tk * (q1*Rd+Rw) )
end





subroutine specific_hum( T, RH, P, q )
! calc. specific humidity : q
! T:  Air temperature ( C. deg.)
! RH:  Relative humidity (%  0-100)
! P:  Air pressure (total)  (hPa)
  implicit none
  real ::  T, RH, P
  real ::  q, r		! mix. ratio
  real ::  e
  real, parameter :: Rd=287.0, Rw = 461.0

  call satvap(T, e)
!   write(*,*) 'e=', e
!  call satvap2(T, e)
!   write(*,*) 'e=', e

  e = (RH/100.0) * e

  r =  e/(P-e) * (Rd/Rw)
  q = r/(1+r)
return
end
  


subroutine satvap(T, e)
! calc. saturation vapor pressure (hPa)  
!  by  Paul.R.Lowe 1977, J. Appl. Met. 16, 100-103
 implicit none
 real :: T   !  temperature in C deg.
 real :: e   !  saturate vapor pressure in hPa

! if (T >= 0.0) then
  e= 6.107799961 + T * (4.436518521e-1 &
                 + T * (1.428945805e-2 &
                 + T * (2.650648471e-4 &
                 + T * (3.031240396e-6 &
                 + T * (2.034080948e-8 &
                 + T * (6.136820929e-11))))))
! else 
!   e= 					&
!   6.109177956 + T * (5.03469897e-1	&
!   + T * (1.886013408e-2		&
!   + T * (4.176223716e-4		&
!   + T * (5.824720280e-6		&
!   + T * (4.838803174e-8		&
!   + T * (1.838826904e-10		&
!   ))))))
!
! end if

return
end

subroutine satvap2(T, e)
! calc. saturation vapor pressure (hPa) by Tetens
 implicit none
 real :: T   !  temperature in C deg.
 real :: e   !  saturate vapor pressure in hPa

! for water
 e = 6.1078*10**(7.5*T/(237.3+T))
! for ice
! e = 6.1078*10**(9.5*T/(265.3+T))

end

subroutine satvap2a(T, e)
! calc. saturation vapor pressure (hPa) by Tetens
 implicit none
 real :: T   !  temperature in C deg.
 real :: e   !  saturate vapor pressure in hPa

! for water
 e = 6.1078*10**(7.5*T/(237.3+T))
! for ice
 e = 6.1078*10**(9.5*T/(265.3+T))

end

subroutine calc_Tdew( e, Tdew )
! calc. dew point temp. from vapor press
 implicit none
 real :: e, Tdew
 real :: f, a, b
   a = 7.5; b=237.3
   f = (log(e)-log(6.1078))/log(10.0)
   Tdew = f * b/(a-f)

  if(Tdew < 0 ) then
   a = 9.5; b=265.3
   f = (log(e)-log(6.1078))/log(10.0)
   Tdew = f * b/(a-f)
  end if 
end

real function dqsatdt( tempC, press )
!   d q_sat/dT
! tempC (C deg.)
!       t (K), press (hPa): Atmospreric pressure

 real :: temp, press, tempC
 real, parameter  ::  eps = 0.622, T_ZERO=273.15
 real :: esat, desdt
 real :: t, ttt

  temp = tempC + T_ZERO
  call satvap2a( tempC, esat)
  t = 1.0 - (1.0 - eps)*(esat/press)
  dqsatdt = ((eps/press)/(t*t))*desdt(temp)
  write(*,*) dqsatdt, ttt

end function

real function desdt( temp )
!  d (esat)/dT
!      esat (hPa)   T (K)
!
 real :: temp ! in K
 real :: tc,  a, b
 real, parameter :: e0=6.1078
 real, parameter  ::  t_ice = 273.15

 tc = temp - t_ice
 if( tc >= 0  ) then
    a = 7.5 ; b=237.3
 else
    a = 9.5; b=265.3
 endif
    desdt = e0 * (10.0**(a*tc/(b+tc))) * (log(10.0) * a*b)/((b+tc)*(b+tc))
 return
end

