! Created by Andrey Debolskiy on 17.12.2023.

module state_utils

    public geo
    public get_geopotential_moist
    public theta2ta
    public ta2theta
    public get_exner
    public get_sigma_from_z
    public get_z_from_sigma
    public get_density_theta
    public get_density
    public dqsdt_sc, qmax_sc, qmax
    public DQSDT2, DQSDT, CLDT
    contains
        subroutine get_geopotential_moist(F,t, qv, sig, fs, rd, eps, kmax)
            implicit none
            integer, intent(in):: kmax
            real, intent(in):: t(kmax)
            real, intent(in):: qv(kmax)
            real, intent(in):: sig(kmax)
            real, intent(in):: rd
            real, intent(in):: eps
            real, intent(in):: fs
            real, intent(out)::f(kmax)

            integer ka, k
            real thalf, qhalf

            F(kmax) = fs - log(sig(kmax)) * rd * t(kmax) * (1.0 + eps * qv(kmax)) / (1.0 + qv(kmax))

            do k = kmax-1, 1, -1
                ka = kmax+1
                thalf = 0.5 * (t(ka) + t(k))
                qhalf = 0.5 * (qv(ka) + qv(k))
                F(k) = F(ka) - (log(sig(ka)) -log(sig(k))) &
                    * rd * thalf * (1.0 + eps * qhalf) / (1.0 + qhalf)
            end do
        end subroutine get_geopotential_moist

        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

    real function get_density( p, t, qv, fluid)
        use phys_fluid, only: fluidParamsDataType
        implicit none
        type(fluidParamsDataType),intent(in):: fluid
        real, intent(in):: p, t, qv

        get_density = (p / (fluid%R_d * t)) * (1.0 + qv) / (1 + qv * fluid%eps_v)
    end function get_density

    real function get_density_theta( p, th, qv, fluid)
        use phys_fluid, only: fluidParamsDataType
        implicit none
        type(fluidParamsDataType),intent(in):: fluid
        real, intent(in):: p, th, qv

        real t
        t = (p / fluid%p0)**(fluid%R_d / fluid%cp)
        get_density_theta = (p / (fluid%R_d * t)) * (1.0 + qv) / (1 + qv * fluid%eps_v)
    end function get_density_theta

    subroutine get_exner( exnr, p, fluid, kmax)
        use phys_fluid, only: fluidParamsDataType
        implicit none
        integer,intent(in):: kmax
        type(fluidParamsDataType), intent(in):: fluid
        real,dimension(kmax), intent(in):: p
        real,dimension(kmax), intent(out):: exnr

        integer k

        do k = kmax,1,-1
            exnr(k) = (p(k)/fluid%p0)**(fluid%R_d / fluid%cp)
        end do
    end subroutine get_exner

    subroutine get_exner_sig( exnr, ps, fluid, grid, kmax)
        use pbl_grid, only : pblgridDataType
        use phys_fluid, only: fluidParamsDataType
        implicit none
        integer,intent(in):: kmax
        type(fluidParamsDataType), intent(in):: fluid
        type(pblgridDataType), intent(in)::grid
        real, intent(in):: ps
        real,dimension(kmax), intent(out):: exnr

        integer k
        do k = kmax,1,-1
            exnr(k) = (ps * grid%sigma(k)/fluid%p0)**(fluid%R_d / fluid%cp)
        end do
    end subroutine get_exner_sig

    subroutine get_sigma_from_z( grid, bl, ps, fluid)
        use pbl_grid, only : pblgridDataType
        use phys_fluid, only: fluidParamsDataType
        use scm_state_data, only: stateBLDataType
        implicit none
        type(fluidParamsDataType), intent(in):: fluid
        type(stateBLDataType), intent(in):: bl
        real, intent(in):: ps
        type(pblgridDataType), intent(inout)::grid


        real dp, rho, p
        integer k, kmax
        kmax = grid%kmax
        grid%sigma_edge(kmax) = 1.0
        rho = get_density(ps, bl%temp(kmax), bl%qv(kmax), fluid)
        dp = - rho * fluid%g * grid%dzc(kmax)
        grid%sigma(kmax) = (ps + dp/2) / ps
        grid%sigma_edge(kmax) = (ps + dp) / ps
        p = ps + dp
        do k = kmax-1,1,-1
            rho = get_density(p, bl%temp(k), bl%qv(k), fluid)
            dp = - rho * fluid%g * grid%dzc(k)
            grid%sigma_edge(k) = grid%sigma_edge(k+1) + dp / ps
            rho = get_density(p, 0.5*(bl%temp(k+1) + bl%temp(k)), &
                    0.5*(bl%qv(k+1) + bl%qv(k)), fluid)
            grid%sigma(k) = grid%sigma(k+1)  &
                    - rho * fluid%g * grid%dze(k)
            p = grid%sigma_edge(k) * ps
        end do
        grid%dsigmac(2:kmax) = grid%sigma(2:kmax) - grid%sigma(1:kmax-1)
        grid%sigma_edge(0) = 0.0
        grid%dsigmae(1:kmax) = grid%sigma_edge(1:kmax) - grid%sigma_edge(1:kmax-1)
    end subroutine get_sigma_from_z

    subroutine get_z_from_sigma( grid, bl, zs, fluid, kmax)
        use pbl_grid, only : pblgridDataType
        use phys_fluid, only: fluidParamsDataType
        use scm_state_data, only: stateBLDataType
        implicit none
        integer,intent(in):: kmax
        type(fluidParamsDataType), intent(in):: fluid
        type(stateBLDataType), intent(in):: bl
        real, intent(in):: zs
        type(pblgridDataType), intent(inout)::grid

        real fs
        real, dimension(kmax):: f
        integer k
        call get_geopotential_moist(f,bl%temp, bl%qv, grid%sigma, &
                fluid%g * zs, fluid%R_d, fluid%eps_v, kmax)
        grid%z_cell(1:kmax) = f(1:kmax) / fluid%g
        grid%z_edge(kmax) = zs
        do k = kmax - 1, 0, -1
            grid%z_edge(k) = grid%z_cell(kmax + 1) + 0.5 *( grid%z_cell(kmax + 1) - grid%z_edge(kmax + 1))
        end do
        grid%dzc(2:kmax) = grid%z_cell(2:kmax) - grid%z_cell(1:kmax-1)
    end subroutine get_z_from_sigma
end module state_utils