model.f90

Go to the documentation of this file.

module model
    use globals, only:dp
    implicit none
    integer :: &
        ngas,&    !< number of gases
        nfv
    integer, dimension(:), allocatable :: &
        mor
    real(dp), dimension(:), allocatable :: &
        dz,&    !< size of finite volume
        dif,&   !< diffusivity
        sol,&   !< solubility
        bc
    private
    public modelpu,ngas,nfv,dz,dif,sol,bc,mor
contains
!********************************BEGINNING*************************************
subroutine modelpu(neq, time, ystate, yprime)
    use constants
    integer, intent(in) :: neq
    real(dp), intent(in) :: time
    real(dp), intent(in) :: ystate(neq)
    real(dp), intent(out) :: yprime(neq)
    integer :: i, j, k
    real(dp) :: fluxw,fluxe
    k=1
    j=1
    do i=1,ngas
        fluxw=-2*dif(k)*sol(k)*(ystate(k)-bc(i))/dz(j)
        fluxe=-2*dif(k)*dif(k+ngas)*sol(k)*sol(k+ngas)*&
            (ystate(k+ngas)-ystate(k))/&
            (dif(k+ngas)*dz(j)*sol(k)+dif(k)*dz(j+1)*sol(k+ngas))
        yprime(k)=(fluxw-fluxe)/dz(j)
        k=k+1
    enddo
    do j=2,nfv-1
        do i=1,ngas
            fluxw=-2*dif(k-ngas)*dif(k)*sol(k)*sol(k-ngas)*&
                (ystate(k)-ystate(k-ngas))/&
                (dif(k)*dz(j-1)*sol(k-ngas)+dif(k-ngas)*dz(j)*sol(k))
            fluxe=-2*dif(k)*dif(k+ngas)*sol(k)*sol(k+ngas)*&
                (ystate(k+ngas)-ystate(k))/&
                (dif(k+ngas)*dz(j)*sol(k)+dif(k)*dz(j+1)*sol(k+ngas))
            yprime(k)=(fluxw-fluxe)/dz(j)
            k=k+1
        enddo
    enddo
    j=nfv
    do i=1,ngas
        fluxw=-2*dif(k-ngas)*dif(k)*sol(k)*sol(k-ngas)*&
            (ystate(k)-ystate(k-ngas))/&
            (dif(k)*dz(j-1)*sol(k-ngas)+dif(k-ngas)*dz(j)*sol(k))
        fluxe=0
        yprime(k)=(fluxw-fluxe)/dz(j)
        k=k+1
    enddo
end subroutine modelpu
!***********************************END****************************************
end module model