mjwl.F File Reference

#include "implicit_f.inc"
#include "mvsiz_p.inc"
#include "com08_c.inc"
#include "param_c.inc"

Go to the source code of this file.

Functions/Subroutines
subroutine	mjwl (pm, mat, off, sig, eint, psh, p0, qold, vol, bfrac, voln, qnew, sold1, sold2, sold3, dvol, nel, df, aburn, er1v, er2v, wdr1v, wdr2v, w1, amu)

Function/Subroutine Documentation

mjwl()

subroutine mjwl	(		pm,
		integer, dimension(mvsiz)	mat,
			off,
			sig,
			eint,
			psh,
			p0,
			qold,
			vol,
			bfrac,
			voln,
			qnew,
			sold1,
			sold2,
			sold3,
			dvol,
		integer	nel,
			df,
			aburn,
			er1v,
			er2v,
			wdr1v,
			wdr2v,
			w1,
			amu )

Definition at line 28 of file mjwl.F.

C-----------------------------------------------
C   I m p l i c i t   T y p e s
C-----------------------------------------------
#include      "implicit_f.inc"
C-----------------------------------------------
C   G l o b a l   P a r a m e t e r s
C-----------------------------------------------
#include      "mvsiz_p.inc"
C-----------------------------------------------
C   D u m m y   A r g u m e n t s
C-----------------------------------------------
      INTEGER MAT(MVSIZ),NEL
      my_real
     .   pm(npropm,*),
     .   off(*),sig(nel,6),eint(*),qold(*),vol(*),bfrac(*),voln(mvsiz),qnew(*),
     .   psh(*),p0(*),
     .   dvol(*),sold1(*),sold2(*),sold3(*), df(*), er1v(*), er2v(*), wdr1v(*), wdr2v(*), w1(*),
     .   aburn(mvsiz),alpha_unit,amu(*)
C-----------------------------------------------
C   C o m m o n   B l o c k s
C-----------------------------------------------
#include      "com08_c.inc"
#include      "param_c.inc"
C-----------------------------------------------
C   L o c a l   V a r i a b l e s
C-----------------------------------------------
      INTEGER I,QOPT
      my_real
     .   einc(mvsiz) , espe(mvsiz), w1df(mvsiz)
      my_real
     .   tbegin, tend,
     .   volo(mvsiz),facm(mvsiz),pold(mvsiz),pnew(mvsiz),
     .   eadd, lambda,rr,rr2,a,m,n,bulk
C-----------------------------------------------
C   S o u r c e   L i n e s
C-----------------------------------------------
      DO i=1,nel
        dvol(i)=half*dvol(i)
        pold(i)=third*(sold1(i)+sold2(i)+sold3(i))
        einc(i)=dvol(i)*(pold(i)-psh(i)-qold(i)-qnew(i))
        eint(i)=eint(i)+einc(i)
        qold(i)=qnew(i)
        volo(i)=voln(i)/df(i)
      ENDDO
 
      !================      AFTERBURNING     ============================!
      ! computes E <- E + lambda(t)*E_add
      qopt       = nint(pm(042,mat(1)))
      eadd       = pm(160,mat(1))
      tbegin     = pm(161,mat(1))
      tend       = pm(162,mat(1))
      rr         = pm(163,mat(1))
      a          = pm(164,mat(1))
      m          = pm(165,mat(1))
      n          = pm(166,mat(1))
      rr2        = pm(167,mat(1))
      alpha_unit = pm(168,mat(1))
      bulk       = pm(044,mat(1))
      IF(eadd==zero)THEN
      !=== no afterburning  
      ! EINT(I) = EINT(I)  
      ELSEIF(qopt==0)THEN
      !=== instantaneous release
        DO i=1,nel
          lambda = zero
          IF(tt > tend .AND. aburn(i)==one)THEN   
            lambda   = one
            aburn(i) = one
          ELSEIF (tt <= tbegin)THEN
            lambda   = zero
            aburn(i) = zero
          ELSE
            lambda   = one
            eint(i)  = eint(i)+(lambda-aburn(i))*eadd*max(em20,volo(i))
            aburn(i) = one
          ENDIF
        ENDDO   
      ELSEIF(qopt==1)THEN
      !=== afterburning with constant rate from Tbegin to Tend
        DO i=1,nel
          lambda = zero
          IF(tt > tend .AND. aburn(i)==one)THEN 
            lambda   = one
            aburn(i) = one
          ELSEIF (tt <= tbegin)THEN
            lambda   = zero
            aburn(i) = zero
          ELSE
            lambda   = (tt-tbegin)*rr
            lambda   = min(one,lambda)
            eint(i)  = eint(i)+(lambda-aburn(i))*eadd*max(em20,volo(i))
            aburn(i) = lambda
          ENDIF          
        ENDDO        
      ELSEIF(qopt==2)THEN
      !=== afterburning with linear rate from Tbegin to Tend
        DO i=1,nel
          lambda = zero
          IF(tt > tend .AND. aburn(i)==one)THEN    ! .AND. ABURN(I)==ONE  needed to add last increment
            lambda   = one
            aburn(i) = one
          ELSEIF (tt <= tbegin)THEN
            lambda   = zero
            aburn(i) = zero
          ELSE
            lambda   = half*rr*tt**2 - rr*tbegin*tt + rr2
            lambda   = max(zero,min(one,lambda))
            eint(i)  = eint(i)+(lambda-aburn(i))*eadd*max(em20,volo(i))
            aburn(i) = lambda
          ENDIF
        ENDDO    
      ELSEIF(qopt==3)THEN
      !=== Miller s extension, rate is depedent on Pressure
        DO i=1,nel
          lambda = zero
          IF(-pold(i)-psh(i) > zero )THEN 
            lambda=aburn(i)+ dt1*a*exp( m*log(one+aburn(i)) )*exp(n*log(alpha_unit*(-pold(i)-psh(i))))
            lambda  = max(lambda,zero)
            lambda  = min(lambda,one)
            eint(i) = eint(i)+(lambda-aburn(i))*eadd*max(em20,volo(i))
            aburn(i)= lambda
          ENDIF
        ENDDO 
      ENDIF
      !================  END AFTERBURNING     ============================!
 
      DO i=1,nel
        espe(i)=eint(i)/max(em20,volo(i))
        w1df(i)=bfrac(i)*w1(i)/df(i)
        facm(i)=bfrac(i)*(wdr1v(i)*er1v(i)+wdr2v(i)*er2v(i))
      ENDDO
 
      DO i=1,nel
        pnew(i)= - psh(i) + (one - bfrac(i)) * (p0(i) + bulk * amu(i)) + 
     .        (facm(i)+(espe(i))*w1df(i))/(one +w1df(i)*dvol(i)/max(em20,volo(i)))
      ENDDO
 
      DO i=1,nel
        !!! Taking -PSH as the minimum value
        pnew(i)= max(zero - psh(i), pnew(i))*off(i)
      ENDDO
 
      DO i=1,nel
        einc(i)= einc(i)-(pnew(i) + psh(i))*dvol(i)
        eint(i)=(eint(i)-(pnew(i) + psh(i))*dvol(i))/max(em20,vol(i))
      ENDDO
 
      DO i=1,nel
        sig(i,1)=sig(i,1)*off(i)-pnew(i)
        sig(i,2)=sig(i,2)*off(i)-pnew(i)
        sig(i,3)=sig(i,3)*off(i)-pnew(i)
      ENDDO
      
C-----------------------------------------------
      RETURN