!
!  $Author: panetta $
!  $Date: 2007/07/20 13:58:43 $
!  $Revision: 1.7 $
!
MODULE Shall_Souza

  USE Constants, ONLY :  &
       cp                 , &
       hl                 , &
       gasr               , &
       grav               , &
       rmwmd              , &
       rmwmdi             , &
       e0c                , &
       delq               , &
       tov                , &
       r8                 , &
       i8
  USE Options, ONLY :       &
       rccmbl            , &
       sthick            , &
       sacum             , &
       acum0             , &
       tbase             , &
       ki                , &
       mlrg              , &
       is                , &
       doprec            , &
       cflric            , &
       ifilt             , &
       dt                , &
       kt                , &
       ktp               , &
       jdt               , &
       nfprt
  IMPLICIT NONE
  PRIVATE
  !!PUBLIC :: InitShall_Souza
  PUBLIC :: shallsouza

  REAL(KIND=r8) :: aa(15)
  REAL(KIND=r8) :: ad(15)
  REAL(KIND=r8) :: ac(15)
  REAL(KIND=r8) :: actop
  REAL(KIND=r8) :: thetae(151,181)
  REAL(KIND=r8) :: tfmthe(431,241)
  REAL(KIND=r8) :: qfmthe(431,241)
  REAL(KIND=r8) :: ess
  INTEGER :: kbase
  INTEGER :: kcr
  REAL(KIND=r8) , ALLOCATABLE :: dels  (:)
  REAL(KIND=r8) , ALLOCATABLE :: gams  (:)
  REAL(KIND=r8) , ALLOCATABLE :: gammod(:)
  REAL(KIND=r8) , ALLOCATABLE :: delmod(:)
  REAL(KIND=r8) :: rlocp
  REAL(KIND=r8) :: rgrav
  REAL(KIND=r8) :: rlrv
  REAL(KIND=r8)  :: const1
  REAL(KIND=r8)  :: const2
  REAL(KIND=r8)  :: xx1

  REAL(KIND=r8) , PARAMETER :: xkapa=0.2857143_r8


CONTAINS


  SUBROUTINE shallsouza(te,qe,ps,sl,tfz,qfz,dt,nCols,kMax,kuo,  &
       noshal, klcl,ktop,par6,par7)
    ! Shallow cumulus heating and moistening tendencies
    ! Enio Pereira de Souza 12/Jul/2001
    ! modified and adapted to GCM by Silvio Nilo Figueroa Jun-2004
    !
    INTEGER, INTENT(IN   ) :: nCols
    INTEGER, INTENT(IN   ) :: kMax
    REAL(KIND=r8) ,    INTENT(INOUT) :: te    (nCols,kMax)  ! air temperature (K)
    REAL(KIND=r8) ,    INTENT(INOUT) :: qe    (nCols,kMax)  ! specific humidity (kg/kg)
    INTEGER, INTENT(INOUT) :: klcl  (nCols)
    INTEGER, INTENT(INOUT) :: ktop  (nCols)
    INTEGER, INTENT(INOUT) :: noshal(nCols)
    REAL(KIND=r8) ,    INTENT(IN   ) :: ps    (nCols)       ! surface pressure (cb)
    REAL(KIND=r8) ,    INTENT(IN   ) :: sl    (kMax)       ! sigma layers
    REAL(KIND=r8) ,    INTENT(IN   ) :: tfz   (nCols)       ! sensible heat flux (W/m^2)
    REAL(KIND=r8) ,    INTENT(IN   ) :: qfz   (nCols)       ! latent heat flux (W/m^2)
    REAL(KIND=r8) ,    INTENT(IN   ) :: dt                 ! time step (s)
    REAL(KIND=r8) ,    INTENT(IN   ) :: par6
    REAL(KIND=r8) ,    INTENT(IN   ) :: par7
    !
    INTEGER :: kzi  (nCols)
    INTEGER :: kuo  (nCols)
    LOGICAL :: llift(nCols)
    LOGICAL :: lconv(nCols)
    !
    REAL(KIND=r8)  :: press   (nCols,kMax)
    REAL(KIND=r8)  :: tv      (nCols,kMax)
    REAL(KIND=r8)  :: ze      (nCols,kMax)
    REAL(KIND=r8)  :: den     (nCols,kMax)
    REAL(KIND=r8)  :: delz    (nCols,kMax)
    REAL(KIND=r8)  :: pi      (nCols,kMax)
    REAL(KIND=r8)  :: th      (nCols,kMax)
    REAL(KIND=r8)  :: thv     (nCols,kMax)
    REAL(KIND=r8)  :: es      (nCols,kMax)
    REAL(KIND=r8)  :: qes     (nCols,kMax)
    REAL(KIND=r8)  :: se      (nCols,kMax)
    REAL(KIND=r8)  :: uhe     (nCols,kMax)
    REAL(KIND=r8)  :: uhes    (nCols,kMax)
    REAL(KIND=r8)  :: uhc     (nCols,kMax)
    REAL(KIND=r8)  :: dtdt    (nCols,kMax)
    REAL(KIND=r8)  :: dqdt    (nCols,kMax)
    REAL(KIND=r8)  :: gamma   (nCols,kMax)
    REAL(KIND=r8)  :: dlamb   (nCols,kMax)
    REAL(KIND=r8)  :: dldzby2 (nCols,kMax)
    REAL(KIND=r8)  :: sc      (nCols,kMax)
    REAL(KIND=r8)  :: sc0     (nCols,kMax)
    REAL(KIND=r8)  :: qc      (nCols,kMax)
    REAL(KIND=r8)  :: ql      (nCols,kMax)
    REAL(KIND=r8)  :: tvm     (nCols,kMax)
    REAL(KIND=r8)  :: scv     (nCols,kMax)
    REAL(KIND=r8)  :: sev     (nCols,kMax)
    REAL(KIND=r8)  :: sc0v    (nCols,kMax)
    REAL(KIND=r8)  :: scvm    (nCols,kMax)
    REAL(KIND=r8)  :: sevm    (nCols,kMax)
    REAL(KIND=r8)  :: sc0vm   (nCols,kMax)
    REAL(KIND=r8)  :: buoy1   (nCols)
    REAL(KIND=r8)  :: buoy2   (nCols)
    REAL(KIND=r8)  :: cape    (nCols)
    REAL(KIND=r8)  :: tcape   (nCols)
    REAL(KIND=r8)  :: emp     (nCols,kMax)
    REAL(KIND=r8)  :: efic    (nCols)
    REAL(KIND=r8)  :: fin     (nCols)
    REAL(KIND=r8)  :: tcold   (nCols)
    REAL(KIND=r8)  :: thot    (nCols)
    REAL(KIND=r8)  :: sigwshb (nCols)
    REAL(KIND=r8)  :: wc      (nCols,kMax)
    REAL(KIND=r8)  :: wssc    (nCols,kMax)
    REAL(KIND=r8)  :: wqsc    (nCols,kMax)
    REAL(KIND=r8)  :: dsdt    (nCols,kMax)
    REAL(KIND=r8)  :: tpar    (nCols)
    REAL(KIND=r8)  :: qex1    (nCols)
    REAL(KIND=r8)  :: qpar    (nCols)
    REAL(KIND=r8)  :: espar   (nCols)
    REAL(KIND=r8)  :: qspar   (nCols)
    REAL(KIND=r8)  :: qexces  (nCols)
    REAL(KIND=r8)  :: dqdp    (nCols)
    REAL(KIND=r8)  :: deltap  (nCols)
    REAL(KIND=r8)  :: plcl    (nCols)
    REAL(KIND=r8)  :: tlcl    (nCols)
    INTEGER    :: icount(nCols)

    !snf
    !REAL :: rho
    !REAL :: dp
    !--------
    !
    ! constants
    !
    !
    REAL(KIND=r8) , PARAMETER :: grav=9.806_r8
    REAL(KIND=r8) , PARAMETER :: vlat=2.5e6_r8
    REAL(KIND=r8) , PARAMETER :: cp=1004.5_r8
    REAL(KIND=r8) , PARAMETER :: rd=287.0_r8
    REAL(KIND=r8) , PARAMETER :: rcp=0.286_r8
    REAL(KIND=r8) , PARAMETER :: p00=100000.0_r8
    REAL(KIND=r8) , PARAMETER :: gammad=0.00976_r8
    REAL(KIND=r8) , PARAMETER :: es00=611.2_r8
    REAL(KIND=r8) , PARAMETER :: epslon=0.622_r8
    REAL(KIND=r8) , PARAMETER :: ummeps=0.378_r8
    REAL(KIND=r8) , PARAMETER :: ta0=273.15_r8
    REAL(KIND=r8) , PARAMETER :: co1=21709759.15_r8
    REAL(KIND=r8) , PARAMETER :: co2=29.65_r8
    REAL(KIND=r8) , PARAMETER :: co3=17.67_r8
    REAL(KIND=r8) , PARAMETER :: c0=0.0_r8
    REAL(KIND=r8) , PARAMETER :: dlamb0=1.0e-6_r8

    REAL(KIND=r8)             :: zref
    REAL(KIND=r8) , PARAMETER :: dthv=2.0_r8
    REAL(KIND=r8) , PARAMETER :: fifty=50.0_r8
    !
    INTEGER    :: i
    INTEGER    :: k
    INTEGER    :: ki

    ! 300 400  560 570 800  was tested
    zref=par6
    !
    DO i=1,nCols
       klcl    (i)=1
       ktop    (i)=1
       noshal  (i)=1
       tcape   (i)=0.0_r8
       sigwshb (i)=0.0_r8
       cape    (i)=0.0_r8
       efic    (i)=0.0_r8
       buoy2   (i)=0.0_r8
       buoy1   (i)=0.0_r8
    END DO

    DO k=1,kMax
       DO i=1,nCols
          dtdt (i,k)=0.0_r8
          dqdt (i,k)=0.0_r8
          dsdt (i,k)=0.0_r8
          ze   (i,k)=0.0_r8
          ql   (i,k)=0.0_r8
          press(i,k)=0.0_r8
          tv   (i,k)=0.0_r8
          th   (i,k)=0.0_r8
       END DO
    END DO
    !
    ! begining of a long loop in index " i "
    !
    !
    ! constructing height profile
    !
    DO k=1,kMax
       DO i=1,nCols
          press(i,k)=ps(i)*sl(k)
          press(i,k)=press(i,k)*1000.0_r8
          tv(i,k)=te(i,k)*(1.0_r8+0.608_r8*qe(i,k))
          th(i,k)=te(i,k)*EXP(rcp*LOG(p00/press(i,k)))
          pi(i,k)=cp*EXP(rcp*LOG(press(i,k)/p00))
          thv(i,k)=th(i,k)*(1.0_r8+0.608_r8*qe(i,k))
          es(i,k)=es00*EXP(co3*(te(i,k)-ta0)/(te(i,k)-co2))
          qes(i,k)=epslon*es(i,k)/(press(i,k)-ummeps*es(i,k))
          !snf added
          IF(qes(i,k) .LE. 1.0e-08_r8) qes(i,k)=1.0e-08_r8
          !-----
          den(i,k)=press(i,k)/(rd*tv(i,k))
       END DO
    END DO

    DO i=1,nCols
       ze(i,1)=29.25_r8*te(i,1)*(101.3_r8-ps(i))/ps(i)
       ze(i,1)=MAX(0.0_r8,ze(i,1))
    END DO

    DO k=2,kMax
       DO i=1,nCols
          delz(i,k)=0.5_r8*rd*(tv(i,k-1)+tv(i,k))* &
               LOG(press(i,k-1)/press(i,k))/grav
          ze(i,k)=ze(i,k-1)+delz(i,k)
       END DO
    END DO
    !!snf added
    !-----------------
    !        DO k=1,kMax
    !        DO i=1,nCols
    !            rho=press(i,k)/(287.0_r8*tv(i,k))
    !            dp=ps(i)*del(k)*1000.0_r8
    !            delz(i,k)=dp/(9.81_r8*rho)
    !        enddo
    !        enddo
    !
    !
    !        DO i=1,nCols
    !            ze(i,1)=delz(i,1)/2.0_r8
    !            ze(i,1)=MAX(0.0_r8,ze(i,1))
    !        dO k=2,kMax
    !            ze(i,k)=ze(i,k-1)+delz(i,k)
    !        END DO
    !        END DO
    !------------
    ki=1
    !
    DO i=1,nCols
       llift(i)=.FALSE.
       qex1(i) = 0.0_r8
       qpar(i) = qe(i,ki)
    END DO
    !
    !
    !     lift parcel from k=ki until it becomes saturated
    !
    DO k = ki, kMax
       DO i = 1, nCols
          IF(kuo(i) /= 1 .or. tfz(i) >= 0.0_r8)THEN
             IF(.NOT.llift(i)) THEN
                tpar(i) = te(i,ki)* &
                     EXP(rcp*LOG(press(i,k)/press(i,ki)))
                espar(i) = es2(tpar(i))
                qspar(i) = epslon*espar(i)*1000.0_r8/ &
                     (press(i,k)-(1-epslon)*espar(i)*1000.0_r8)
                qexces(i) = qpar(i) - qspar(i)
                !
                !     if parcel not saturated,  try next level
                !
                IF (qexces(i).LT.0.0_r8) THEN
                   qex1(i) = qexces(i)
                   !
                   !     saturated - go down and find p,t, sl at the lcl;
                   !     if sat. exists in first layer (k=ki), use this as lcl
                   !
                ELSE IF (k.EQ.ki) THEN
                   plcl(i)  = press(i,k)
                   tlcl(i)  = tpar(i)
                   tlcl(i)  = tlcl(i) + 1.0_r8
                   klcl(i)  = k
                   llift(i) = .TRUE.
                ELSE
                   dqdp(i)   = (qexces(i)-qex1(i))/ &
                        (press(i,k-1)-press(i,k))
                   deltap(i) = qexces(i)/dqdp(i)
                   plcl(i)   = press(i,k) + deltap(i)
                   tlcl(i)   = tpar(i) * (1.0_r8+2.0_r8*rcp*deltap(i) &
                        /(press(i,k)+press(i,k-1)))
                   tlcl(i)   = tlcl(i) + 1.0_r8
                   klcl(i)   = k
                   llift(i) = .TRUE.
                   !
                   !     give parcel a one-degree goose
                   !
                END IF
                !
                !     lifting cond level found - get out of loop
                !
             END IF
          END IF
       END DO
    END DO
    !
    ! testing the difference in height between lcl and zi, if zi < z(lcl)
    ! there will be no shallow cumulus
    !
    kzi=0
    DO k=3,kMax
       DO i=1,nCols
          IF((thv(i,k)-thv(i,k-1)) > dthv .AND. kzi(i) == 0 )THEN
             kzi(i)=k
          END IF
       END DO
    END DO
    DO i=1,nCols
       lconv(i) = kuo(i) == 1 .or. tfz(i) <= 0.0_r8 .or. .NOT.llift(i) .or. &
            kzi(i) < klcl(i) .or. klcl(i) <= 2
    END DO
    !
    !     quit if parcel still unsat at k = kthick - - -
    !     in this case,set low value of plcl as signal to shalmano
    !
    DO k=1,kMax
       DO i=1,nCols
          IF (.NOT.lconv(i)) THEN
             !
             ! static energy profiles
             !
             se   (i,k)=cp*te(i,k)+grav*ze(i,k)
             uhe  (i,k)=se(i,k)+vlat*qe(i,k)
             uhes (i,k)=se(i,k)+vlat*qes(i,k)
             gamma(i,k)=co1*press(i,k)*qes(i,k)*qes(i,k)/es(i,k)
             gamma(i,k)=gamma(i,k)/(te(i,k)*te(i,k))
          END IF
       END DO
    END DO

    DO k=2,MAXVAL(klcl)
       DO i=1,nCols
          IF (.NOT.lconv(i) .AND. k <= klcl(i) ) THEN
             !
             ! vertical profile of the entrainment rate and cloud moist
             ! static energy
             !
             uhc(i,1)=uhe(i,1)
             uhc(i,k)=uhe(i,1)
          END IF
       END DO
    END DO

    DO k=MINVAL(klcl),kMax
       DO i=1,nCols
          IF (.NOT.lconv(i).AND. k >= klcl(i)+1) THEN
             dlamb  (i,k)=EXP(LOG(dlamb0)+2.3_r8*ze(i,k)/zref)
             dlamb  (i,k)=MIN(0.1_r8,dlamb(i,k))
             dldzby2(i,k)=dlamb(i,k)*delz(i,k)/2.0_r8
             uhc    (i,k)=(uhc(i,k-1)-dldzby2(i,k)*(uhc(i,k-1)-uhe(i,k)- &
                  uhe(i,k-1)))/(1+dldzby2(i,k))
          END IF
       END DO
    END DO

    !
    ! calculating cloud variables qc, ql, sc
    !
    DO k=1,kMax
       DO i=1,nCols
          IF (.NOT.lconv(i)) THEN
             sc (i,k)=se (i,k)+(uhc(i,k)-uhes(i,k))/(1+gamma(i,k))
             sc0(i,k)=se (i,k)+(uhc(i,1)-uhes(i,k))/(1+gamma(i,k))
             qc (i,k)=qes(i,k)+gamma(i,k)*(uhc(i,k)-uhes(i,k))/ &
                  (vlat*(1+gamma(i,k)))
             ql (i,k)=0.0_r8
          END IF
       END DO
    END DO

    DO k=MINVAL(klcl),kMax
       DO i=1,nCols
          IF (.NOT.lconv(i) .AND. k >= klcl(i)+1) THEN
             ql(i,k)= ql(i,k-1)-(qc(i,k)-qc(i,k-1))-dlamb(i,k)* &
                  (qc(i,k)-qe(i,k))*delz(i,k)- &
                  (c0+dlamb(i,k))*ql(i,k-1)*delz(i,k)
             ql(i,k)=MAX(0.00000001_r8,ql(i,k))
          END IF
       END DO
    END DO

    !
    ! determining cloud top based on integrated buoyancy
    !
    DO k=1,kMax
       DO i=1,nCols
          IF (.NOT.lconv(i)) THEN
             scv (i,k)=sc (i,k)+cp*te(i,k)*(0.608_r8*qc(i,k)-ql(i,k))
             sev (i,k)=se (i,k)+0.608_r8*cp*te(i,k)*qe(i,k)
             sc0v(i,k)=sc0(i,k)+0.608_r8*cp*te(i,k)*qe(i,k)
          END IF
       END DO
    END DO

    DO k=2,kMax
       DO i=1,nCols
          IF (.NOT.lconv(i)) THEN
             scvm (i,k)=(scv(i,k)+scv(i,k-1))/2.0_r8
             sevm (i,k)=(sev(i,k)+sev(i,k-1))/2.0_r8
             sc0vm(i,k)=(sc0v(i,k)+sc0v(i,k-1))/2.0_r8
             tvm  (i,k)=(tv(i,k)+tv(i,k-1))/2.0_r8
          END IF
       END DO
    END DO

    DO i=1,nCols
       IF (.NOT.lconv(i)) THEN
          !
          ! determination of the integrated buoyancy between surface and lcl
          ! The calculation assumes that the surface flux is in W/m2. Therefore
          ! we divide it by density*cp in order to convert it to Km/s
          !
          buoy1(i)=0.0_r8
          buoy1(i)=tfz(i)*(1.0_r8+0.608_r8*qe(i,1))/(cp*den(i,1))
          buoy1(i)=grav*ze(i,klcl(i))*buoy1(i)/tv(i,1)
          buoy1(i)=EXP(0.29_r8+0.6667_r8*LOG(buoy1(i)))
          !
          ! checking wether the parcel is able to sustain positive
          ! buoyancy one level above lcl
          !
          cape(i)=0.0_r8
          buoy2(i)=0.0_r8
          buoy2(i)=buoy2(i)+gammad*(scvm(i,klcl(i)+1)-sevm(i,klcl(i)+1))* &
               delz(i,klcl(i)+1)/tvm(i,klcl(i)+1)
       END IF
    END DO

    DO i=1,nCols
       lconv(i) = lconv(i) .or. (buoy1(i)+buoy2(i)) <= 0.0_r8
    END DO



    DO k=MINVAL(klcl),kMax
       DO i=1,nCols
          IF (.NOT.lconv(i) .and. k >= klcl(i)+2) THEN
             !
             ! calculating cloud top and cape
             !
             buoy2(i)=buoy2(i)+gammad*(scvm(i,k)-sevm(i,k))* &
                  delz(i,k)/tvm(i,k)

             IF ((buoy1(i) + buoy2(i)) <= 0.0_r8 .AND. ktop(i) == 1) THEN
                ktop(i)=k-1
             END IF

          END IF
       END DO
    END DO
    DO k=MINVAL(klcl),MAXVAL(ktop)
       DO i=1,nCols
          IF (.NOT.lconv(i) .AND. k >= klcl(i) .AND. k <= ktop(i)) THEN
             emp (i,k)=(sc0vm(i,k)-sevm(i,k))
             emp (i,k)=MAX(0.0_r8,emp(i,k))
             cape(i  )=cape(i)+gammad*emp(i,k)*delz(i,k)/tvm(i,k)
          END IF
       END DO
    END DO
    !
    ! calculating the cloud base mass flux
    !
    ! se quer aumentar fluxo de massa usar tcold(i)=te(i,klcl)
    !
    DO i=1,nCols
       IF (.NOT.lconv(i)) THEN
          thot (i) = te(i,1)
          tcold(i) = te(i,klcl(i))
       END IF
    END DO

    !    icount=1
    !    DO i=1,nCols
    !      IF (.NOT.lconv(i)) THEN
    !        DO k=3,ktop(i)
    !          tcold (i)=tcold (i)+te(i,k)
    !          icount(i)=icount(i)+1
    !        END DO
    !      END IF
    !    END DO

    icount=1
    DO k=3,MAXVAL(ktop)
       DO i=1,nCols
          IF (.NOT.lconv(i) .AND.k <= ktop(i)) THEN
             tcold (i)=tcold (i)+te(i,k)
             icount(i)=icount(i)+1
          END IF
       END DO
    END DO

    DO i=1,nCols
       IF (.NOT.lconv(i)) THEN
          tcold(i)=tcold(i)/icount(i)
          efic(i)=(thot(i)-tcold(i))/thot(i)
       END IF
    END DO

    DO i=1,nCols
       lconv(i)=lconv(i).or.efic(i).LE.0.0_r8.OR.cape(i).LT.40.0_r8
    END DO

    DO i=1,nCols
       IF (.NOT.lconv(i)) THEN
          !
          fin(i)=tfz(i)+qfz(i)
          !
          !snf
          !       tcape parece esta forte......  trocar...
          !       tcape(i)=cape(i)
          !snf    tcape(i)=2*cape(i)
          !       tcape(i)=1.25*cape(i)
          !GOOD   tcape(i)=1.6 *cape(i)
          !
          tcape(i)=par7*cape(i)
          !
          ! noshall parameter define where shallow convection works
          !
          noshal(i)=0
          !
          ! fluxo na base da nuvem
          !
          sigwshb(i)=efic(i)*fin(i)/(den(i,klcl(i))*tcape(i))
       END IF
    END DO

    wssc=0.0_r8
    wqsc=0.0_r8

    DO k= MINVAL(klcl),MAXVAL(ktop)
       DO i=1,nCols
          IF (.NOT.lconv(i) .AND. k >= klcl(i) .AND. k <= ktop(i)) THEN
             wc(i,k)=sigwshb(i)*((ze(i,ktop(i))-ze(i,k))/ &
                  (ze(i,ktop(i))-ze(i,klcl(i))))
             wssc(i,k)=wc(i,k)*(sc(i,k)-vlat*ql(i,k)-se(i,k))
             wqsc(i,k)=wc(i,k)*(qc(i,k)+ql(i,k)-qe(i,k))
          END IF
       END DO
    END DO

    DO k=MINVAL(klcl),MAXVAL(ktop)
       DO i=1,nCols
          IF (.NOT.lconv(i) .AND. k >= klcl(i)+1 .AND. k <= ktop(i)-1) THEN
             dsdt(i,k)=-(wssc(i,k+1)-wssc(i,k-1))/(ze(i,k+1)-ze(i,k-1))
             dtdt(i,k)=  dsdt(i,k)/pi(i,k)
             dqdt(i,k)=-(wqsc(i,k+1)-wqsc(i,k-1))/(ze(i,k+1)-ze(i,k-1))
          END IF
       END DO
    END DO
    !
    ! updating temperature and moisture fields due to shallow convection
    !
    DO k=1,kMax
       DO i=1,nCols
          te(i,k)=te(i,k)+dtdt(i,k)*2.0_r8*dt
          qe(i,k)=qe(i,k)+dqdt(i,k)*2.0_r8*dt
       END DO
    END DO

  END SUBROUTINE shallsouza

  REAL(KIND=r8)  FUNCTION es2(t)
    REAL(KIND=r8) , INTENT(IN) :: t
    REAL(KIND=r8)  :: tx
    REAL(KIND=r8) , PARAMETER  :: d1=0.6107042e0_r8
    REAL(KIND=r8) , PARAMETER  :: d2=4.441157e-2_r8
    REAL(KIND=r8) , PARAMETER  :: d3=1.432098e-3_r8
    REAL(KIND=r8) , PARAMETER  :: d4=2.651396e-5_r8
    REAL(KIND=r8) , PARAMETER  :: d5=3.009998e-7_r8
    REAL(KIND=r8) , PARAMETER  :: d6=2.008880e-9_r8
    REAL(KIND=r8) , PARAMETER  :: d7=6.192623e-12_r8

    ! statement function

    tx = t - tbase
    IF (tx >= -50.0_r8) THEN
       es2 = d1 + tx*(d2 + tx*(d3 + tx*(d4 + tx*(d5 + tx*(d6 + d7*tx)))))
    ELSE
       es2=0.00636e0_r8*EXP(25.6e0_r8*(tx+50.e0_r8)/(tbase-50.e0_r8))
    END IF
  END FUNCTION es2

END MODULE Shall_Souza
