Sources Esope | Cast3M

mshear
C MSHEAR    SOURCE    PV        11/03/07    21:17:33     6885
      SUBROUTINE MSHEAR(DDEP,FOR0,FORF,VAR0,VARF,NVAR,DDINL,
     & DTRANP,DTRANM,BETA,NPELA,TRFA,DOCP,DOCM,EXPN,
     & CURFP,NCURFP,CURKP,NCURKP,CURLP,NCURLP,
     & CURFM,NCURFM,CURKM,NCURKM,CURLM,NCURLM,KERRE)
      IMPLICIT INTEGER(I-N)
      IMPLICIT REAL*8(A-H,O-Z)
C
C======================================================================
C       MUR_SHEAR, Elisa, Armelle, Alessandra et Pierre, ELSA/ISPRA
C                  01/2002 --> ../2004 ...
C======================================================================
C
C                    MODELE GLOBAL DE MUR EN CISAILLEMENT
C           (Sur des elements de poutre TIMO - Effort tranchant/Cisail.)
C
C======================================================================
C=======================================================================
C
C LISTE D'ECHANGE
C ---------------
C
C DDEP         = Incrément de déplacement
C FOR0         = Effort initial
C FORF         = Effort final
C VAR0         = Variables internes initiales
C VARF         = Variables internes finales
C NVAR         = Nb de variable interne
C DDINL        = Increment de deformation inelastique
C
C DTRANP       = zone de transition  (sens positif)
C DTRANM       =                     (sens négatif)
C
C BETA         = Coefficient BETA (equivalence nrj-deplacement)
C NPELA        = Nb de pts par unite de zone elastique pour les
C                segments non lineaires
C+2003
C TRFA         = Coefficient d'entrainement des origine (>0,&lt;1)
C DOCP         = Deplacement d'ouverture des fissure + (>0)
C DOCM         = Deplacement d'ouverture des fissure - (&lt;0)
C+2003
C+2004
C EXPN         = Exposant n de la fonction de cumul du domage
C+2004
C
C CURFP,NCURFP = courbe de charge  +           (x>0,y>0)
C CURKP,NCURKP = courbe de raideur +           (x>0,y>0)
C CURLP,NCURLP = courbe de domaine elastique + (x>0,y>0)
C CURFM,NCURFM = courbe de charge  -           (x&lt;0,y&lt;0)
C CURKM,NCURKM = courbe de raideur -           (x&lt;0,y>0)
C CURLM,NCURLM = courbe de domaine elastique - (x&lt;0,y>0)
C
C VARIABLES INTERNES (0=initiales, F=finales)
C ------------------
C
C        : Numero du cas
C DEP    : Deplacement total
C EPLUS  : Energie "positive"
C EMOIN  : Energie "negative"
C DPLUS  : Deplacement de reference +
C DMOIN  : Deplacement de reference -
C FCINI  : Force de reference pour la surface cinematique
C FCAMP  : Taille de la surface cinematique
C KCINE  : Raideur de la surface cinematique
C APLUS  : Coefficient multiplicateur de l'amplitude +
C AMOIN  : Coefficient multiplicateur de l'amplitude -
C+2003
C OPLUS  : Position de l'origine + &lt;0
C OMOIN  : Position de l'origine - &lt;0
C+2003
C
C
C=======================================================================
C
      PARAMETER (XZER=0.D0,UN=1.D0,XPETIT=1.D-7)
      DIMENSION CURFP(2,NCURFP),CURKP(2,NCURKP),CURLP(2,NCURLP)
      DIMENSION CURFM(2,NCURFM),CURKM(2,NCURKM),CURLM(2,NCURLM)
C
      REAL*8 KCINE0,KCINEF
C
      DIMENSION VAR0(NVAR),VARF(NVAR)
      LOGICAL LSWICH
C
      KERRE=0
C
C  Lecture variables internes
C
      ICAS  = nint(VAR0(1))
      DEP0  = VAR0(2)
      EPLUS0= VAR0(3)
      EMOIN0= VAR0(4)
      DPLUS0= VAR0(5)
      DMOIN0= VAR0(6)
      FCINI0= VAR0(7)
      FCAMP0= VAR0(8)
      KCINE0= VAR0(9)
      APLUS0= VAR0(10)
      AMOIN0= VAR0(11)
C+2003
      OPLUS0= VAR0(12)
      OMOIN0= VAR0(13)
C+2003
C+2004
      E= VAR0(12)
      OMOIN0= VAR0(13)
C+2004
C
C  Initialisation des variables finales
C
      DDINL=XZER
      DEPF  =DEP0+DDEP
      DO IE1=1,NVAR
        VARF(IE1)=VAR0(IE1)
      ENDDO
      VARF(2)=DEPF
      FORF=FOR0
*
*     WRITE(6,*)ICAS,DEP0,DEPF
*
C
C  Cas de l'increment nul
C
      IF(DDEP.EQ.XZER)RETURN
C==================================================================
C
C                          DRIVER
C
C==================================================================
      GOTO(1,2,3,4,5,6,7),ICAS+1
C==========================================================
C   A: CAS DU MODELE DANS SON ETAT ELASTIQUE        ICAS= 0
C==========================================================
 1    CONTINUE
C-->Condition de branchement sur B0 ou C0
      IF(DEPF.GT.CURFP(1,2))THEN
        ICAS=1
        DEP0   = CURFP(1,2)
        FOR0   = CURFP(2,2)
        KCINE0 = FOR0/DEP0
        DMOIN0 = CURFM(1,2)
        VARF(6)= DMOIN0
        GOTO 2
      ELSEIF(DEPF.LT.CURFM(1,2))THEN
        ICAS=4
        DEP0   = CURFM(1,2)
        FOR0   = CURFM(2,2)
        KCINE0 = FOR0/DEP0
        DPLUS0 = CURFP(1,2)
        VARF(5)= DPLUS0
        GOTO 5
      ENDIF
C-->Regime normal: On calcul l'etat elastique
      IF(DEPF.GE.0)THEN
        KCINEF=CURFP(2,2)/CURFP(1,2)
        FORF=KCINEF*DEPF
      ELSE
        KCINEF=CURFM(2,2)/CURFM(1,2)
        FORF=KCINEF*DEPF
      ENDIF
      RETURN
C==========================================================
C   B0: CAS DU CHARGEMENT MONOTONE DIRECTION +      ICAS= 1
C==========================================================
 2    CONTINUE
C-->Condition de branchement sur C1
      IF(DDEP.LT.XZER)THEN
        ICAS=2
        GOTO 3
      ENDIF
C-->Regime normal: On procede sur la courbe de charge +
C 2003DPLUSF=DEPF
      DPLUSF=DEPF-OPLUS0
C 2003CALL YOFXCU(DEPF  ,CURFP,NCURFP,  FORF,DYSDX,KERRE)
      CALL YOFXCU(DPLUSF,CURFP,NCURFP,  FORF,DYSDX,KERRE)
      FORF=FORF*(UN-APLUS0)
      VARF(1)=ICAS
      VARF(5)=DPLUSF
C-->On calcule et on sauve les caracteristiques de la surface
C   Cinematique en tenant compte de la zone de transition
C 2003CALL MSHEA1(DEPF,ICAS,FORF,DPLUSF,DMOIN0, DTRANP,DTRANM,
      CALL MSHEA1(DPLUSF,ICAS,FORF,DPLUSF,DMOIN0, DTRANP,DTRANM,
     >            APLUS0,CURFP,NCURFP,CURKP,NCURKP,CURLP,NCURLP,
     >            AMOIN0,CURFM,NCURFM,CURKM,NCURKM,CURLM,NCURLM,
     >            FCINIF,FCAMPF,KCINEF)
C2003>            FCINIF,FCAMPF,KCINEF,ALPHP)
      VARF(7) =FCINIF
      VARF(8) =FCAMPF
      VARF(9) =KCINEF
C-->Increment de d.p.
      DDINL=DDINL+DEPF-FORF/KCINEF-(DEP0-FOR0/KCINE0)
C+2003
      IF(DPLUSF.GE.DOCP)THEN
        IF(OMOIN0.EQ.XZER)THEN
          OMOINF=TRFA*(DPLUSF-DOCP)
        ELSE
          OMOINF=OMOIN0+TRFA*(DPLUSF-DPLUS0)
        ENDIF
        VARF(13)=OMOINF
      ENDIF
C+2003
      RETURN
C==========================================================
C   C1: BRANCHE DE UNLOADING ELASTIQUE -            ICAS= 2
C==========================================================
 3    CONTINUE
C-->Condition de branchement sur B0, B2 ou C2
      FORF=FOR0+KCINE0*DDEP
      IF(FORF.GT.FCINI0)THEN
C -->On sort par le haut (B0 ou B2)
        DEP0=DEP0+(FCINI0-FOR0)/KCINE0
        DDEP=DEPF-DEP0
        FOR0=FCINI0
C  -->On retourne sur B0... si on vient de B0
        IF(ABS(DEP0-DPLUS0).LT.XPETIT*DPLUS0)THEN
          ICAS=1
          GOTO 2
C  -->On retourne sur B2... si on ne vient PAS de B0
        ELSE
          ICAS=6
          GOTO 7
        ENDIF
      ELSEIF(FORF.LT.FCINI0-FCAMP0)THEN
C --> On continue sur C2
        DEP0=DEP0+(FCINI0-FCAMP0-FOR0)/KCINE0
        DDEP=DEPF-DEP0
        FOR0=FCINI0-FCAMP0
C
C+2002 On repointe "tangent" a la courbe (si possible)
        CALL YOFXCU(DMOIN0,CURFM,NCURFM,  FMOIN0,DYSDX,KERRE)
C 2003  STRIAL=(FMOIN0-FOR0)/(DMOIN0-DEP0)
        STRIAL=(FMOIN0-FOR0)/(DMOIN0+OMOIN0-DEP0)
        DO IE1=1,NCURFM
          IF(DMOIN0.GT.CURFM(1,IE1))THEN
C 2003      SNEXT=(CURFM(2,IE1)-FOR0)/(CURFM(1,IE1)-DEP0)
            SNEXT=(CURFM(2,IE1)-FOR0)/(CURFM(1,IE1)+OMOIN0-DEP0)
            IF(SNEXT.GT.STRIAL)THEN
              DMOIN0=CURFM(1,IE1)
              STRIAL=SNEXT
            ELSE
              GOTO 31
            ENDIF
          ENDIF
        ENDDO
 31     CONTINUE
        VARF(6)=DMOIN0
C+2002
        ICAS=3
        GOTO 4
      ENDIF
C-->Regime normal: On decharge elastiquement
      VARF(1)=ICAS
      RETURN
C==========================================================
C   C2: BRANCHE DE UNLOADING ANELASTIQUE -          ICAS= 3
C==========================================================
 4    CONTINUE
C-->Condition de branchement sur B1
      IF(DDEP.GT.XZER)THEN
        ICAS=5
        GOTO 6
      ENDIF
C-->Valeur de ALPHP
C 2003ALPHP=(DEP0-FOR0/KCINE0-DTRANM)/(DTRANP-DTRANM)
C 2003ALPHP=MAX(XZER,ALPHP)
C 2003ALPHP=MIN(UN  ,ALPHP)
C+2003
CC    CALL MSHEA2(DEP0,DTRANP+OPLUS0,DTRANM+OMOIN0,
      CALL MSHEA2(DEP0,DTRANP+OPLUS0+OMOIN0,DTRANM+OPLUS0+OMOIN0,
     >            FOR0,FCAMP0,KCINE0,ICAS,
     >                 ALPHP,ALPHM)
C+2003
C-->Determination de la taille du sous increment
C   WARNING DDEP et XINCR sont negatifs
      XINCR=(CURFM(2,2)-CURFP(2,2))/NPELA
      NINCR=INT(DDEP/XINCR)+1
      XINCR=DDEP/NINCR
C-->Determination initiale du point "parfait" que l'on vise (DMOIN0,FMOIN0)
      CALL YOFXCU(DMOIN0,CURFM,NCURFM,  FMOIN0,DYSDX,KERRE)
C-->Loop sur les sous-increments
      LSWICH=.FALSE.
      DO IE1=1,NINCR
C-->Traitement du branchement possible sur C0
C 2003  IF(DEP0+XINCR.LT.DMOIN0)THEN
C       IF(DEP0+XINCR.LT.DMOIN0+OMOIN0)THEN
C 2003  IF(DEP0+XINCR.LT.(UN-XPETIT)*DMOIN0)THEN
        IF(DEP0+XINCR.LT.(UN-XPETIT)*(DMOIN0+OMOIN0))THEN
          LSWICH=.TRUE.
C 2003    XINCR =DMOIN0-DEP0
          XINCR =DMOIN0+OMOIN0-DEP0
        ENDIF
C-->On incremente explicitement tout ce qui doit etre incremente
C   WARNING:force effective pointee=FMOIN0*(UN-AMOIN0)
C 2003  FORF=FOR0+XINCR/(DMOIN0-DEP0)*(FMOIN0*(UN-AMOIN0)-FOR0)
        FORF=FOR0+XINCR/(DMOIN0+OMOIN0-DEP0)*(FMOIN0*(UN-AMOIN0)-FOR0)
        DEP0=DEP0+XINCR
        DENE=ABS(FCAMP0/2*(XINCR-(FORF-FOR0)/KCINE0))
        DDINL=DDINL+XINCR+FOR0/KCINE0
        FOR0=FORF
C
        EPLUS0=EPLUS0+   ALPHP  *DENE
C 2003  EMOIN0=EMOIN0+(UN-ALPHP)*DENE
        EMOIN0=EMOIN0+   ALPHM  *DENE
*       APLUS0=MIN(UN,APLUS0+BETA*   ALPHP*  DENE)
*       AMOIN0=MIN(UN,AMOIN0+BETA*(UN-ALPHP)*DENE)
C 2003  APLUS0=TANH(BETA*EPLUS0)
C 2003  AMOIN0=TANH(BETA*EMOIN0)
        CALL MSHEA3(BETA,EXPN, EPLUS0,EMOIN0, APLUS0,AMOIN0)
C-->On calcule les caracteristiques de la surface
C   Cinematique en tenant compte de la zone de transition
        CALL MSHEA1(DEP0,ICAS,FOR0,DPLUS0,DMOIN0,
     >              DTRANP,DTRANM,
     >              APLUS0,CURFP,NCURFP,CURKP,NCURKP,CURLP,NCURLP,
     >              AMOIN0,CURFM,NCURFM,CURKM,NCURKM,CURLM,NCURLM,
     >              FCINI0,FCAMP0,KCINE0)
C2003>              FCINI0,FCAMP0,KCINE0,ALPHP)
C+2003
CC      CALL MSHEA2(DEP0,DTRANP+OPLUS0,DTRANM+OMOIN0,
        CALL MSHEA2(DEP0,DTRANP+OPLUS0+OMOIN0,DTRANM+OPLUS0+OMOIN0,
     >              FOR0,FCAMP0,KCINE0,ICAS,
     >                   ALPHP,ALPHM)
C+2003
        DDINL=DDINL-FOR0/KCINE0
C-->Switch sur C0
        IF(LSWICH)THEN
          ICAS=4
          VARF(3)=EPLUS0
          VARF(4)=EMOIN0
          VARF(10)=APLUS0
          VARF(11)=AMOIN0
          GOTO 5
        ENDIF
C-->Fin normale loop sur les sous-increments
      ENDDO
C
      FORF=FOR0
C
      VARF(1)=ICAS
      VARF(3)=EPLUS0
      VARF(4)=EMOIN0
      VARF(7)=FCINI0
      VARF(8)=FCAMP0
      VARF(9)=KCINE0
      VARF(10)=APLUS0
      VARF(11)=AMOIN0
C
      RETURN
C==========================================================
C   B0: CAS DU CHARGEMENT MONOTONE DIRECTION -      ICAS= 4
C==========================================================
 5    CONTINUE
C-->Condition de branchement sur B1
      IF(DDEP.GT.XZER)THEN
        ICAS=5
        GOTO 6
      ENDIF
C-->Regime normal: On procede sur la courbe de charge -
C 2003DMOINF=DEPF
      DMOINF=DEPF-OMOIN0
C 2003CALL YOFXCU(DEPF  ,CURFM,NCURFM,  FORF,DYSDX,KERRE)
      CALL YOFXCU(DMOINF,CURFM,NCURFM,  FORF,DYSDX,KERRE)
      FORF=FORF*(UN-AMOIN0)
      VARF(1)=ICAS
      VARF(6)=DMOINF
C-->On calcule et on sauve les caracteristiques de la surface
C   Cinematique en tenant compte de la zone de transition
C 2003CALL MSHEA1(DEPF,ICAS,FORF,DPLUS0,DMOINF, DTRANP,DTRANM,
      CALL MSHEA1(DEPF,ICAS,FORF,DPLUS0,DMOINF, DTRANP,DTRANM,
     >            APLUS0,CURFP,NCURFP,CURKP,NCURKP,CURLP,NCURLP,
     >            AMOIN0,CURFM,NCURFM,CURKM,NCURKM,CURLM,NCURLM,
     >            FCINIF,FCAMPF,KCINEF)
C2003>            FCINIF,FCAMPF,KCINEF,ALPHP)
      VARF(7) =FCINIF
      VARF(8) =FCAMPF
      VARF(9) =KCINEF
C-->Increment de d.p.
      DDINL=DDINL+DEPF-FORF/KCINEF-(DEP0-FOR0/KCINE0)
C+2003
      IF(DMOINF.LE.DOCM)THEN
        IF(OPLUS0.EQ.XZER)THEN
          OPLUSF=TRFA*(DMOINF-DOCM)
        ELSE
          OPLUSF=OPLUS0+TRFA*(DMOINF-DMOIN0)
        ENDIF
        VARF(12)=OPLUSF
      ENDIF
C+2003
      RETURN
C==========================================================
C   B1: BRANCHE DE RELOADING ELASTIQUE +            ICAS= 5
C==========================================================
 6    CONTINUE
C-->Condition de branchement sur C0, C2 ou B2
      FORF=FOR0+KCINE0*DDEP
      IF(FORF.LT.FCINI0)THEN
C -->On sort par le bas (C0 ou C2)
        DEP0=DEP0+(FCINI0-FOR0)/KCINE0
        DDEP=DEPF-DEP0
        FOR0=FCINI0
C  -->On retourne sur C0... si on vient de C0
        IF(ABS(DEP0-DMOIN0).LT.ABS(XPETIT*DMOIN0))THEN
          ICAS=4
          GOTO 5
C  -->On retourne sur C2... si on ne vient PAS de C0
        ELSE
          ICAS=3
          GOTO 4
        ENDIF
      ELSEIF(FORF.GT.FCINI0+FCAMP0)THEN
C --> On continue sur B2
        DEP0=DEP0+(FCINI0+FCAMP0-FOR0)/KCINE0
        DDEP=DEPF-DEP0
        FOR0=FCINI0+FCAMP0
C
C2002 On repointe "tangent" a la courbe (si possible)
        CALL YOFXCU(DPLUS0,CURFP,NCURFP,  FPLUS0,DYSDX,KERRE)
C 2003  STRIAL=(FPLUS0-FOR0)/(DPLUS0-DEP0)
        STRIAL=(FPLUS0-FOR0)/(DPLUS0+OPLUS0-DEP0)
        DO IE1=1,NCURFP
          IF(DPLUS0.LT.CURFP(1,IE1))THEN
C 2003      SNEXT=(CURFP(2,IE1)-FOR0)/(CURFP(1,IE1)-DEP0)
            SNEXT=(CURFP(2,IE1)-FOR0)/(CURFP(1,IE1)+OPLUS0-DEP0)
            IF(SNEXT.GT.STRIAL)THEN
              DPLUS0=CURFP(1,IE1)
              STRIAL=SNEXT
            ELSE
              GOTO 61
            ENDIF
          ENDIF
        ENDDO
 61     CONTINUE
        VARF(5)=DPLUS0
C2002
        ICAS=6
        GOTO 7
      ENDIF
C-->Regime normal: On decharge elastiquement
      VARF(1)=ICAS
      RETURN
C==========================================================
C   B2: BRANCHE DE UNLOADING ANELASTIQUE -          ICAS= 6
C==========================================================
 7    CONTINUE
C-->Condition de branchement sur C1
      IF(DDEP.LT.XZER)THEN
        ICAS=2
        GOTO 3
      ENDIF
C-->Valeur de ALPHP
C 2003ALPHP=(DEP0-FOR0/KCINE0-DTRANM)/(DTRANP-DTRANM)
C 2003ALPHP=MAX(XZER,ALPHP)
C 2003ALPHP=MIN(UN  ,ALPHP)
C+2003
CC    CALL MSHEA2(DEP0,DTRANP+OPLUS0,DTRANM+OMOIN0,
      CALL MSHEA2(DEP0,DTRANP+OPLUS0+OMOIN0,DTRANM+OPLUS0+OMOIN0,
     >            FOR0,FCAMP0,KCINE0,ICAS,
     >                 ALPHP,ALPHM)
C+2003
 
C-->Determination de la taille du sous increment
C   WARNING DDEP et XINCR sont negatifs
      XINCR=(CURFP(2,2)-CURFM(2,2))/NPELA
      NINCR=INT(DDEP/XINCR)+1
      XINCR=DDEP/NINCR
C-->Determination initiale du point "parfait" que l'on vise (DPLUS0,FPLUS0)
      CALL YOFXCU(DPLUS0,CURFP,NCURFP,  FPLUS0,DYSDX,KERRE)
C-->Loop sur les sous-increments
      LSWICH=.FALSE.
      DO IE1=1,NINCR
C-->Traitement du branchement possible sur B0
C 2003  IF(DEP0+XINCR.GT.DPLUS0)THEN
C       IF(DEP0+XINCR.GT.DPLUS0+OPLUS0)THEN
C 2003  IF(DEP0+XINCR.GT.(UN-XPETIT)*DPLUS0)THEN
        IF(DEP0+XINCR.GT.(UN-XPETIT)*(DPLUS0+OPLUS0))THEN
          LSWICH=.TRUE.
C 2003    XINCR =DPLUS0-DEP0
          XINCR =DPLUS0+OPLUS0-DEP0
        ENDIF
C-->On incremente explicitement tout ce qui doit etre incremente
C   WARNING:l'increment de d.p. est calcule en 2 steps
C   WARNING:force effective pointee=FPLUS0*(UN-APLUS0)
C 2003  FORF=FOR0+XINCR/(DPLUS0-DEP0)*(FPLUS0*(UN-APLUS0)-FOR0)
        FORF=FOR0+XINCR/(DPLUS0+OPLUS0-DEP0)*(FPLUS0*(UN-APLUS0)-FOR0)
        DEP0=DEP0+XINCR
        DENE=ABS(FCAMP0/2*(XINCR-(FORF-FOR0)/KCINE0))
        DDINL=DDINL+XINCR+FOR0/KCINE0
        FOR0=FORF
C
        EPLUS0=EPLUS0+   ALPHP  *DENE
C 2003  EMOIN0=EMOIN0+(UN-ALPHP)*DENE
        EMOIN0=EMOIN0+   ALPHM  *DENE
*       APLUS0=MIN(UN,APLUS0+BETA*   ALPHP*  DENE)
*       AMOIN0=MIN(UN,AMOIN0+BETA*(UN-ALPHP)*DENE)
C 2003  APLUS0=TANH(BETA*EPLUS0)
C 2003  AMOIN0=TANH(BETA*EMOIN0)
        CALL MSHEA3(BETA,EXPN, EPLUS0,EMOIN0, APLUS0,AMOIN0)
C-->On calcule les caracteristiques de la surface
C   Cinematique en tenant compte de la zone de transition
        CALL MSHEA1(DEP0,ICAS,FOR0,DPLUS0,DMOIN0,
     >              DTRANP,DTRANM,
     >              APLUS0,CURFP,NCURFP,CURKP,NCURKP,CURLP,NCURLP,
     >              AMOIN0,CURFM,NCURFM,CURKM,NCURKM,CURLM,NCURLM,
     >              FCINI0,FCAMP0,KCINE0)
C2003>              FCINI0,FCAMP0,KCINE0,ALPHP)
C+2003
CC      CALL MSHEA2(DEP0,DTRANP+OPLUS0,DTRANM+OMOIN0,
        CALL MSHEA2(DEP0,DTRANP+OPLUS0+OMOIN0,DTRANM+OPLUS0+OMOIN0,
     >              FOR0,FCAMP0,KCINE0,ICAS,
     >                   ALPHP,ALPHM)
C+2003
        DDINL=DDINL-FOR0/KCINE0
C-->Switch sur B0
        IF(LSWICH)THEN
          ICAS=1
          VARF(3)=EPLUS0
          VARF(4)=EMOIN0
          VARF(6)=DMOIN0
          VARF(10)=APLUS0
          VARF(11)=AMOIN0
          GOTO 2
        ENDIF
C-->Fin normale loop sur les sous-increments
      ENDDO
C
      FORF=FOR0
C
      VARF(1)=ICAS
      VARF(3)=EPLUS0
      VARF(4)=EMOIN0
      VARF(7)=FCINI0
      VARF(8)=FCAMP0
      VARF(9)=KCINE0
      VARF(10)=APLUS0
      VARF(11)=AMOIN0
C
      RETURN
C
      END