! Created by Andrey Debolskiy on 17.12.2023.

module state_utils

    public geo
    public theta2ta
    public ta2theta
    public dqsdt_sc, qmax_sc, qmax
    public DQSDT2, DQSDT, CLDT
    contains

        SUBROUTINE GEO(T,FS,F,SIG, R, LL)
        implicit none
        integer, intent(in):: LL
        real, intent(in):: T(LL)
        real, intent(in):: SIG(LL)
        real, intent(in):: FS
        real, intent(out):: F(LL)
        real, intent(in):: R

        integer k, ka

        F(LL)=FS-ALOG(SIG(LL))*R*T(LL)

        do  K=LL-1,1,-1

            KA=K+1
            F(K)=F(KA)+0.5E0*R*(ALOG(SIG(KA))-ALOG(SIG(K)))* &
             (T(KA)+T(K))

        end do

        END SUBROUTINE GEO




        subroutine theta2ta(ta, theta, p0, sig, appa, kl)
            implicit none
            real, dimension(kl), intent(in):: theta
            real, dimension(kl), intent(in)::sig
            real, intent(in):: p0, appa
            integer, intent(in):: kl
            real, dimension(kl), intent(out):: ta

            integer k

            do k=1,kl
                ta(k) = theta(k) * (p0 * sig(k))**appa
            end do
        end subroutine theta2ta

        subroutine ta2theta(theta, ta, p0, sig, appa, kl)
            implicit none
            real, dimension(kl), intent(in):: ta
            real, dimension(kl), intent(in)::sig
            real, intent(in):: p0, appa
            integer, intent(in):: kl
            real, dimension(kl), intent(out):: theta

            integer k

            do k=1,kl
                theta(k) = ta(k) * (p0 * sig(k))**(-appa)
            end do
        end subroutine ta2theta

        real function qmax_sc(t, p, aa)
            implicit none

            real, intent(in) :: t, p, aa
            real :: pb2, pb

            pb2 = t - 273.2e0
            if ((aa.gt.1.0e0.and.pb2.ge.0.0e0).or.(aa.lt.1.0e0.and.pb2.ge.-12.0e0)) then
                pb=17.55e0/(t-31.1e0)
            else
                pb=21.85e0/(t-7.65e0)
            end if
            qmax_sc=3.80042e-3*exp(pb*pb2)/p
        end function qmax_sc
    real function qmax(T,P, AA )
        IMPLICIT NONE
        real, intent(in):: t, p, aa
        !local
        real pb2, pb
        PB2=T-273.2E0
        IF((AA.GT.1.0E0.AND.PB2.GE.0.0E0).OR. &
                 (AA.LT.1.0E0.AND.PB2.GE.-12.0E0))THEN
            PB=17.55E0/(T-31.1E0 )
        ELSE
            PB=21.85E0/(T-7.65E0)
        ENDIF
        QMAX=3.80042E-3*EXP(PB*PB2)/P
        RETURN
    end function qmax


    real function dqsdt_sc(a, b, h)
            real, intent(in) ::  a, b, h

            dqsdt_sc=4248.855e0*h*qmax_sc(a,b, 0.0e0)/(a-31.10e0)**2
        end function dqsdt_sc

    REAL FUNCTION CLDT(Q,T,P,C)
        !
        IMPLICIT NONE
        REAL, INTENT(IN) ::  Q,T,P,C
        !      REAL DQSDT
        !REAL :: QMAX

        CLDT=(Q-C*QMAX(T,P, 0.0E0))/(1.0E0+2488.851E0*DQSDT(T,P,C))
        RETURN
    END FUNCTION CLDT

    REAL FUNCTION DQSDT(A,B,H)
        !
        IMPLICIT NONE
        REAL, INTENT(IN) ::  A, B, H
        !REAL :: QMAX

        DQSDT=4248.855E0*H*QMAX(A,B, 0.0E0)/(A-31.10E0)**2
        RETURN
    END FUNCTION DQSDT

    REAL FUNCTION DQSDT2(A,B,H)

        IMPLICIT NONE

        REAL, INTENT(IN) ::  A, B, H
        !      REAL  DQSDT
        !REAL :: QMAX

        DQSDT2=4248.855E0*H*(DQSDT(A,B,    H )/(A-31.10E0)**2 - &
                2.0E0*QMAX (A,B, 0.0E0)/(A-31.10E0)**3)
        RETURN
    END FUNCTION DQSDT2

end module state_utils