Skip to content
Snippets Groups Projects
state_utils.f90 3.4 KiB
Newer Older
  • Learn to ignore specific revisions
  • Debolskiy Andrey's avatar
    Debolskiy Andrey committed
    ! 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
    
    Debolskiy Andrey's avatar
    Debolskiy Andrey committed
        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
    
    
    Debolskiy Andrey's avatar
    Debolskiy Andrey committed
    end module state_utils