qsat_util.F90

Go to the documentation of this file.

MODULE qsat_util

!This module contains the subroutines used to generate the lookup table for computing saturation vapour pressure.
!Variables designated as real*8 should be kept at least as so, even if ESTBLX is required in real*4.

IMPLICIT NONE

PRIVATE
PUBLIC esinit

!Saturation vapor pressure table ESTBLX initialization.
!To be defined for entire module:
integer, parameter :: degsubs    =  100
real(8), parameter :: tmintbl    =  150.0, tmaxtbl = 333.0
integer, parameter :: tablesize  =  nint(tmaxtbl-tmintbl)*degsubs + 1

CONTAINS

subroutine esinit(ESTBLX)

 IMPLICIT NONE

 !OUTPUT
 real(8), dimension(TABLESIZE) :: estblx

 !LOCALS
 real(8), parameter :: zeroc = 273.16, tmix = -20.0

 real(8), dimension(TABLESIZE) :: estble, estblw

 integer :: i
 real(8)    :: t, delta_t

 delta_t = 1.0/degsubs

 do i=1,tablesize

    t = (i-1)*delta_t + tmintbl

    if(t>zeroc) then
       call qsatlqu0(estble(i),t)
    else
       call qsatice0(estble(i),t)
    end if

    call qsatlqu0(estblw(i),t)

    t = t-zeroc
    if(t>=tmix .and. t<0.0) then
       estblx(i) = ( t/tmix )*( estble(i) - estblw(i) ) + estblw(i)
    else
       estblx(i) = estble(i)
    end if

 end do

 end subroutine esinit


subroutine qsatlqu0(QS,TL)
!SUPERSATURATED AS LIQUID

 IMPLICIT NONE

 !INPUTS
 real(8) :: TL!, TMAXTBL

 !OUTPUTS
 real(8) :: QS

 !LOCALS
 real(8), parameter :: ZEROC   = 273.16
 real(8), parameter :: TMINLQU = zeroc - 40.0

 real*8,  parameter :: B6 = 6.136820929e-11*100.0
 real*8,  parameter :: B5 = 2.034080948e-8 *100.0
 real*8,  parameter :: B4 = 3.031240396e-6 *100.0
 real*8,  parameter :: B3 = 2.650648471e-4 *100.0
 real*8,  parameter :: B2 = 1.428945805e-2 *100.0
 real*8,  parameter :: B1 = 4.436518521e-1 *100.0
 real*8,  parameter :: B0 = 6.107799961e+0 *100.0

 real(8) :: TX, EX, TI, TT

 tx = tl

 if    (tx<tminlqu) then
    ti = tminlqu
 elseif(tx>tmaxtbl) then
    ti = tmaxtbl
 else
    ti = tx
 end if

 tt = ti-zeroc  !Starr polynomial fit
 ex = (tt*(tt*(tt*(tt*(tt*(tt*b6+b5)+b4)+b3)+b2)+b1)+b0)

 tl = tx
 qs = ex

 return

end subroutine qsatlqu0


subroutine qsatice0(QS,TL)
!SUPERSATURATED AS ICE

 IMPLICIT NONE   
      
 !INPUTS
 real(8) :: TL

 !OUTPUTS
 real(8) :: QS

 !LOCALS
 real(8), parameter :: ZEROC = 273.16, tminstr = -95.0
 real(8), parameter :: TMINICE = zeroc + tminstr

 real(8), parameter :: TSTARR1 = -75.0, tstarr2 = -65.0, tstarr3 = -50.0,  tstarr4 = -40.0

 real*8,  parameter :: BI6= 1.838826904e-10*100.0
 real*8,  parameter :: BI5= 4.838803174e-8 *100.0
 real*8,  parameter :: BI4= 5.824720280e-6 *100.0
 real*8,  parameter :: BI3= 4.176223716e-4 *100.0
 real*8,  parameter :: BI2= 1.886013408e-2 *100.0
 real*8,  parameter :: BI1= 5.034698970e-1 *100.0
 real*8,  parameter :: BI0= 6.109177956e+0 *100.0
 real*8,  parameter :: S16= 0.516000335e-11*100.0
 real*8,  parameter :: S15= 0.276961083e-8 *100.0
 real*8,  parameter :: S14= 0.623439266e-6 *100.0
 real*8,  parameter :: S13= 0.754129933e-4 *100.0
 real*8,  parameter :: S12= 0.517609116e-2 *100.0
 real*8,  parameter :: S11= 0.191372282e+0 *100.0
 real*8,  parameter :: S10= 0.298152339e+1 *100.0
 real*8,  parameter :: S26= 0.314296723e-10*100.0
 real*8,  parameter :: S25= 0.132243858e-7 *100.0
 real*8,  parameter :: S24= 0.236279781e-5 *100.0
 real*8,  parameter :: S23= 0.230325039e-3 *100.0
 real*8,  parameter :: S22= 0.129690326e-1 *100.0
 real*8,  parameter :: S21= 0.401390832e+0 *100.0
 real*8,  parameter :: S20= 0.535098336e+1 *100.0

 real(8) :: TX, TI, TT, W, EX

 tx = tl

 if (tx<tminice) then
    ti = tminice
 elseif(tx>zeroc  ) then
    ti = zeroc
 else
    ti = tx
 end if

 tt = ti - zeroc
 if (tt < tstarr1 ) then
     ex = (tt*(tt*(tt*(tt*(tt*(tt*s16+s15)+s14)+s13)+s12)+s11)+s10)
 elseif(tt >= tstarr1 .and. tt < tstarr2) then
     w = (tstarr2 - tt)/(tstarr2-tstarr1)
     ex =       w *(tt*(tt*(tt*(tt*(tt*(tt*s16+s15)+s14)+s13)+s12)+s11)+s10) &
              + (1.-w)*(tt*(tt*(tt*(tt*(tt*(tt*s26+s25)+s24)+s23)+s22)+s21)+s20)
 elseif(tt >= tstarr2 .and. tt < tstarr3) then
     ex = (tt*(tt*(tt*(tt*(tt*(tt*s26+s25)+s24)+s23)+s22)+s21)+s20)
 elseif(tt >= tstarr3 .and. tt < tstarr4) then
     w = (tstarr4 - tt)/(tstarr4-tstarr3)
     ex =       w *(tt*(tt*(tt*(tt*(tt*(tt*s26+s25)+s24)+s23)+s22)+s21)+s20) &
              + (1.-w)*(tt*(tt*(tt*(tt*(tt*(tt*bi6+bi5)+bi4)+bi3)+bi2)+bi1)+bi0)
 else
     ex = (tt*(tt*(tt*(tt*(tt*(tt*bi6+bi5)+bi4)+bi3)+bi2)+bi1)+bi0)
 endif

 qs = ex

 return
 
end subroutine qsatice0

end MODULE qsat_util