Télécharger la plaquette Cast3M


fr en

Télécharger propha.eso

Retour à la liste

Numérotation des lignes : 

propha

  1. C PROPHA    SOURCE    CB215821  24/04/12    21:16:58     11897          
  2.  
  3.       SUBROUTINE PROPHA
  4.  
  5.       IMPLICIT INTEGER(I-N)
  6.       IMPLICIT REAL*8(A-H,O-Z)
  7. *
  8. *  Calcule la nouvelle proportion de phase
  9. * en entree
  10. *           : son chamelem de materiau associé
  11. *           : un champ de temperature proposée avec les mult de
  12. *             lagrange
  13. *  en sortie MCHAML de proportion de phase mis a jour
  14. *
  15.  
  16.  
  17. -INC PPARAM
  18. -INC CCOPTIO
  19. -INC CCREEL
  20.  
  21. -INC SMMODEL
  22. -INC SMCHAML
  23.       POINTEUR mchelp.mchelm
  24.       POINTEUR mchamp.mchaml
  25.       POINTEUR melvap.melval
  26. -INC SMCHPOI
  27. -INC SMELEME
  28. -INC SMCOORD
  29.  
  30.       SEGMENT SCPR
  31.         INTEGER            ICPR1(NODE)
  32.         INTEGER            ICPR2(NODE)
  33.       ENDSEGMENT
  34.  
  35.       SEGMENT              ICPR3(NODE,N_CON)
  36.       SEGMENT              ICPR4(NODE,N_CON)
  37.  
  38.       SEGMENT S_CONS
  39.         CHARACTER*(LCONMO) CHACON(N_CON)
  40.         INTEGER            LIS_LX(N_CON)
  41.         INTEGER            NB_MOD(N_CON)
  42.         INTEGER            LISMOD(N_CON,N_CON)
  43. C       N_CON : Dimensionnement au plus grand du SEGMENT
  44. C       CHACON: Nom des constituants
  45. C       LIS_LX: Liste de  pointeurs S_LX
  46. C       NB_MOD: Nombre de IMODEL partageant le meme CONSTITUANT
  47. C       LISMOD: Liste des IMODEL partageant le meme CONSTITUANT
  48.       ENDSEGMENT
  49.  
  50.       SEGMENT S_LX
  51.         INTEGER NB_LX
  52.         INTEGER LISELE(N_LX)
  53.         INTEGER LISNOL(N_LX)
  54.         INTEGER LISNOI(N_LX)
  55.         REAL*8  Q_REA (N_LX)
  56.         REAL*8  Q_MAX (N_LX)
  57. C       N_LX  : Dimensionnement au plus grand du SEGMENT
  58. C       NB_LX : Nombre de noeuds 'LX' du CHAMP qui concerne ce constituant
  59. C       LISELE: Liste de Pointeur S_ELEM
  60. C       LISNOL: Numeros des noeuds 'LX'
  61. C       LISNOI: Numeros des noeuds 'INCO' relies aux 'LX'
  62. C       Q_REA : Chaleur de la reaction a distribuer   sur les noeuds
  63. C       Q_MAX : Chaleur maximale         distribuable sur les noeuds
  64.       ENDSEGMENT
  65.  
  66.       SEGMENT S_ELEM
  67.         INTEGER NB_EL
  68.         INTEGER MELPRO(N_EL)
  69.         INTEGER LISNEL(N_EL)
  70.         INTEGER LISNNO(N_EL)
  71.         REAL*8  QTOTNO(N_EL)
  72.         REAL*8  QDEJNO(N_EL)
  73.         REAL*8  QRESNO(N_EL)
  74.         REAL*8  QTOTEL(N_EL)
  75.         REAL*8  QDEJEL(N_EL)
  76.         REAL*8  QRESEL(N_EL)
  77. C       NB_EL : Nombre d'elements de ce type touchant le noeud 'INCO'
  78. C       N_EL  : Dimensionnement au plus grand du SEGMENT
  79. C       MELPRO: Liste des pointeurs MELVAL associes a la proportion de phase
  80. C       LISNEL: Liste des numeros d'elements dans le MELEME SIMPLE
  81. C       LISNNO: Liste des numeros de noeud dans l'element trouve
  82. C       QTOTNO: Chaleur totale         integree au noeud 'INCO'
  83. C       QDEJNO: Chaleur deja accumulee integree au noeud 'INCO'
  84. C       QRESNO: Chaleur restante       integree au noeud 'INCO'
  85. C       QTOTEL: Chaleur totale         sommee dans l'element
  86. C       QDEJEL: Chaleur deja accumulee sommee dans l'element
  87. C       QRESEL: Chaleur restante       sommee dans l'element
  88.       ENDSEGMENT
  89.  
  90.       SEGMENT ITRAV
  91.         INTEGER IPERM1(NNN)
  92.         INTEGER IPERM2(NNN)
  93.         INTEGER IT1(NNN)
  94.         REAL*8  TT1(NNN)
  95.         REAL*8  TT2(NNN)
  96.       ENDSEGMENT
  97.  
  98.       SEGMENT ISEGSU(0)
  99.       SEGMENT ISEG(0)
  100.  
  101.       CHARACTER*(LOCOMP) MOCOMP
  102.       LOGICAL     LOG1
  103.  
  104. C Un peu de marge sur la precision des calculs (pb aix)
  105.       XZPREL=XZPREC*1.D3
  106. * +--------------------------------------------------------------------+
  107. * |    LECTURE DES ARGUMENTS GIBIANE                                   |
  108. * +--------------------------------------------------------------------+
  109.       ipmodl = 0
  110.       NNN    = 0
  111.  
  112. C     Lecture du modele
  113.       call LIROBJ('MMODEL  ',mmodel,1,iretou)
  114.       CALL ACTOBJ('MMODEL  ',mmodel,1)
  115.       if(ierr.ne.0) return
  116.  
  117. *     Lecture du mchaml de proportions de phase (aux NOEUDS)
  118.       call LIROBJ('MCHAML  ',IPIN,1,iretou)
  119.       call ACTOBJ('MCHAML  ',IPIN,1)
  120.       IF(IERR .NE. 0) RETURN
  121.  
  122.       CALL REDUAF(IPIN,mmodel,mchel1,0,IR,KER)
  123.       IF(IR   .NE. 1) CALL ERREUR(KER)
  124.       IF(IERR .NE. 0) RETURN
  125.  
  126.       call quesup(ipmodl,mchel1,0,0,isupp,iret)
  127.       if (isupp.gt.1) then
  128.         call erreur(609)
  129.         return
  130.       endif
  131.  
  132. *     Lecture du mchaml de chaleur latente integree (aux NOEUDS)
  133.       call LIROBJ('MCHAML  ',IPIN,1,iretou)
  134.       call ACTOBJ('MCHAML  ',IPIN,1)
  135.       IF(IERR .NE. 0) RETURN
  136.       CALL REDUAF(IPIN,mmodel,MCHEL2,0,IR,KER)
  137.       IF(IR   .NE. 1) CALL ERREUR(KER)
  138.       IF(IERR .NE. 0) RETURN
  139.       call quesup(ipmodl,MCHEL2,0,0,isupq,iret)
  140.       if (isupq.gt.1) then
  141.         call erreur(609)
  142.         return
  143.       endif
  144.  
  145. C     Lecture du CHPOINT des temperatures proposees
  146.       CALL LIROBJ('CHPOINT ',mchpoi,1,IRETOU)
  147.       CALL ACTOBJ('CHPOINT ',mchpoi,1)
  148.       IF(IERR.NE.0) RETURN
  149.  
  150. C     Initialisation du chapeau du resultat
  151.       SEGINI,MCHELP=MCHEL1
  152.  
  153.       XUN = 1.D0
  154.  
  155. * +--------------------------------------------------------------------+
  156. * |    Reperage du MSOUPO des 'LX' dans mchpoi                         |
  157. * +--------------------------------------------------------------------+
  158.       NSOUPO=mchpoi.ipchp(/1)
  159.       DO iii=1,NSOUPO
  160.         MSOUPO=mchpoi.ipchp(iii)
  161.         NC    =MSOUPO.NOCOMP(/2)
  162.         DO ICOMP=1,NC
  163.           IF(MSOUPO.NOCOMP(ICOMP) .EQ. 'LX  ') THEN
  164.             MPOVAL = IPOVAL
  165.             MELEME = IGEOC
  166.             N      = MPOVAL.VPOCHA(/1)
  167.             GOTO 10
  168.           ENDIF
  169.         ENDDO
  170.       ENDDO
  171.  
  172. C     Pas de 'LX' trouve on sort
  173.       CALL ACTOBJ('MCHAML  ',MCHELP,1)
  174.       CALL ECROBJ('MCHAML  ',MCHELP)
  175.       RETURN
  176.  
  177.  10   CONTINUE
  178.  
  179. C     Si 'LX' trouve, il est la seule composante de son MSOUPO
  180.       IF (ICOMP .NE. 1) CALL ERREUR(5)
  181.  
  182. * +----------------------------------------------------------------------------+
  183. * |    Liste des CONSTITUANTS et 'LX' pour les modeles 'CHANGEMENT_PHASE'      |
  184. * +----------------------------------------------------------------------------+
  185.       NODE  = nbpts
  186.       SEGINI,ISEGSU
  187.       SEGINI,SCPR
  188.       ISEGSU(**)= SCPR
  189.       N1        = MMODEL.KMODEL(/1)
  190.       N_CON     = N1
  191.       SEGINI,S_CONS
  192.       ISEGSU(**)=S_CONS
  193.       SEGINI,ICPR3
  194.       ISEGSU(**)=ICPR3
  195.  
  196.       NB_CON = 0
  197.       DO 1 I_MODE=1,N1
  198. C       Boucle sur les SOUS-MODELES 'CHANGEMENT_PHASE'
  199.         IMODEL=MMODEL.KMODEL(I_MODE)
  200.         NFOR  =IMODEL.FORMOD(/2)
  201.         CALL PLACE(IMODEL.FORMOD,NFOR,IPLAC,'CHANGEMENT_PHASE')
  202.         IF  (IPLAC .EQ. 0) GOTO 1
  203.  
  204. C       Recherche du noeud 'LX' dans le MAILLAGE de BLOCAGE
  205.         IF(TYMODE(2) .NE. 'MAILLAGE')THEN
  206.           CALL ERREUR(5)
  207.           RETURN
  208.         ENDIF
  209.  
  210.         IPT1  = IMODEL.IVAMOD(2)
  211.         NB_EL1= IPT1.NUM(/2)
  212. C       Listes des CONSTITUANTS
  213.         DO 2 I_CONS=1,NB_CON
  214.           IF(S_CONS.CHACON(I_CONS) .EQ. IMODEL.CONMOD)THEN
  215. C           Constituant existant
  216.             nbmode=S_CONS.NB_MOD(I_CONS)+1
  217.             S_CONS.NB_MOD(I_CONS)=nbmode
  218.             S_CONS.LISMOD(nbmode,I_CONS)=IMODEL
  219.             S_LX = S_CONS.LIS_LX(I_CONS)
  220.             N_LX = S_LX.LISELE(/1) + NB_EL1
  221.             SEGADJ,S_LX
  222.             GOTO 3
  223.           ENDIF
  224.  2      CONTINUE
  225. C       Nouvelle Zone trouvee
  226.         NB_CON = NB_CON + 1
  227.         I_CONS = NB_CON
  228.         N_LX   = NB_EL1
  229.         SEGINI,S_LX
  230.         ISEGSU(**)=S_LX
  231.         S_LX.NB_LX = 0
  232.         S_CONS.LIS_LX(I_CONS)  = S_LX
  233.         S_CONS.CHACON(I_CONS)  = IMODEL.CONMOD
  234.         S_CONS.NB_MOD(I_CONS)  = 1
  235.         S_CONS.LISMOD(1,I_CONS)= IMODEL
  236.  3      CONTINUE
  237.  
  238.         DO 4 I_EL1=1,NB_EL1
  239. C         Noeud 1 : 'LX'
  240. C         Noeud 2 : 'Inconnues classiques'
  241.           I_LX   = IPT1.NUM(1,I_EL1)
  242.           I_INCO = IPT1.NUM(2,I_EL1)
  243.           IF(ICPR3(I_INCO,I_CONS) .EQ. 0)THEN
  244.             ICPR3(I_INCO,I_CONS) = I_LX
  245.           ENDIF
  246.           SCPR.ICPR1(I_LX) = I_INCO
  247.           SCPR.ICPR2(I_LX) = I_CONS
  248.  4      CONTINUE
  249.  1    CONTINUE
  250.  
  251.       IF (NB_CON .EQ. 0) THEN
  252. C       Cas ou aucune SOUS-ZONE ne correspond aux modeles 'CHANGEMENT_PHASE'
  253.         CALL ACTOBJ('MCHAML  ',MCHELP,1)
  254.         CALL ECROBJ('MCHAML  ',MCHELP)
  255.         RETURN
  256.       ENDIF
  257.  
  258. * +--------------------------------------------------------------------+
  259. * |    Creation du resultat (extension du MELVAL obligatoire)          |
  260. * +--------------------------------------------------------------------+
  261.       N1=MCHELP.ICHAML(/1)
  262.       DO iN1=1,N1
  263.         MCHAML=MCHELP.ICHAML(iN1)
  264.         IPT2  =MCHELP.IMACHE(iN1)
  265.         N2    =1
  266.         SEGINI,MCHAMP
  267.         MCHELP.ICHAML(iN1)= MCHAMP
  268.         MCHAMP.NOMCHE(1)  ='PPHA'
  269.         MCHAMP.TYPCHE(1)  ='REAL*8'
  270.  
  271.         N1PTEL=IPT2.NUM(/1)
  272.         N1EL  =IPT2.NUM(/2)
  273.         N2PTEL=0
  274.         N2EL  =0
  275.         SEGINI,MELVAP
  276.         MCHAMP.IELVAL(1)=MELVAP
  277.  
  278.         MELVA1=MCHAML.IELVAL(1)
  279.         MAXP1=MELVA1.VELCHE(/1)
  280.         MAXE1=MELVA1.VELCHE(/2)
  281.         DO IE=1,N1EL
  282.           DO IP=1,N1PTEL
  283.             MELVAP.VELCHE(IP,IE)=MELVA1.VELCHE(MIN(IP,MAXP1),
  284.      &                                         MIN(IE,MAXE1))
  285.           ENDDO
  286.         ENDDO
  287.       ENDDO
  288.  
  289. * +--------------------------------------------------------------------+
  290. * |    Reperage des Elements qui touchent un noeud 'INCO'              |
  291. * +--------------------------------------------------------------------+
  292.       SEGINI,ICPR4
  293.       ISEGSU(**)=ICPR4
  294.       DO I_CONS=1,NB_CON
  295.         nbmode=S_CONS.NB_MOD(I_CONS)
  296.         DO I_MODE=1,nbmode
  297.           IMODEL= S_CONS.LISMOD(I_MODE,I_CONS)
  298.           IPT1  = IMODEL.IMAMOD
  299.           DO I_EL1=1,IPT1.NUM(/2)
  300.             DO I_NO1=1,IPT1.NUM(/1)
  301.               NU_INC=IPT1.NUM(I_NO1,I_EL1)
  302.               ICPR4(NU_INC,I_CONS)=ICPR4(NU_INC,I_CONS) + 1
  303.             ENDDO
  304.           ENDDO
  305.  
  306.           IPT2  = IMODEL.IVAMOD(2)
  307.           DO I_EL2=1,IPT2.NUM(/2)
  308.  
  309.           ENDDO
  310.         ENDDO
  311.  
  312. C       Creation des SEGMENTS S_ELEM et placement dans ICPR4
  313.         DO 20 iii=1,NODE
  314.           N_EL=ICPR4(iii,I_CONS)
  315.           IF(N_EL .EQ. 0) GOTO 20
  316.           NNN =MAX(NNN,N_EL)
  317.           SEGINI,S_ELEM
  318.           ISEGSU(**)=S_ELEM
  319.           ICPR4(iii,I_CONS)=S_ELEM
  320.           S_ELEM.NB_EL = 0
  321.  20     CONTINUE
  322.  
  323.         DO I_MODE=1,nbmode
  324.           IMODEL = S_CONS.LISMOD(I_MODE,I_CONS)
  325.           IPT1   = IMODEL.IMAMOD
  326.           NB_NO1 = IPT1.NUM(/1)
  327.           NB_EL1 = IPT1.NUM(/2)
  328.  
  329. C         Recuperation des MELVAL de 'PPHA' et inconnue duale 'Q' de cette SOUS-ZONE
  330.           NOMID  = IMODEL.LNOMID(2)
  331.           MOCOMP  = NOMID.LESOBL(1)
  332.  
  333. *         +--------------------------------------------------------+
  334. *         |    Recherche du MELVAL correspondant a 'PPHA'          |
  335. *         +--------------------------------------------------------+
  336.           DO 41 I_MCH=1,MCHELP.ICHAML(/1)
  337.             IF(MCHELP.CONCHE(I_MCH) .NE. IMODEL.CONMOD) GOTO 41
  338.             IF(MCHELP.IMACHE(I_MCH) .EQ. IPT1  ) THEN
  339.               MCHAM1=MCHELP.ICHAML(I_MCH)
  340.               GOTO 42
  341.             ENDIF
  342.    41     CONTINUE
  343.           CALL ERREUR(472)
  344.           RETURN
  345.    42     CONTINUE
  346.  
  347.           DO 46 I_COMP=1,MCHAM1.NOMCHE(/2)
  348.             IF(MCHAM1.NOMCHE(I_COMP).EQ.'PPHA') THEN
  349.               IF (MCHAM1.TYPCHE(I_COMP).NE.'REAL*8') THEN
  350.                 MOTERR(1:16)  = MCHAM1.TYPCHE(I_COMP)
  351.                 MOTERR(17:20) ='PPHA'
  352.                 MOTERR(21:36) = MCHELP.TITCHE
  353.                 CALL ERREUR(552)
  354.                 RETURN
  355.               ENDIF
  356.               MELVA1=MCHAM1.IELVAL(I_COMP)
  357.               NPT1=MELVA1.VELCHE(/1)
  358.               NEL1=MELVA1.VELCHE(/2)
  359.               GOTO 47
  360.             ENDIF
  361.    46     CONTINUE
  362.           MOTERR = 'PPHA'
  363.           CALL ERREUR(677)
  364.           RETURN
  365.    47     CONTINUE
  366.  
  367. *         +-------------------------------------------------------------+
  368. *         |    Recherche du MELVAL correspondant aux quantites integrees|
  369. *         +-------------------------------------------------------------+
  370.           DO 411 I_MCH=1,MCHEL2.ICHAML(/1)
  371.             IF(MCHEL2.CONCHE(I_MCH) .NE. IMODEL.CONMOD) GOTO 411
  372.             IF(MCHEL2.IMACHE(I_MCH) .EQ. IPT1  ) THEN
  373.               MCHAM2=MCHEL2.ICHAML(I_MCH)
  374.               GOTO 421
  375.             ENDIF
  376.   411     CONTINUE
  377.           CALL ERREUR(472)
  378.           RETURN
  379.   421     CONTINUE
  380.  
  381.           DO 461 I_COMP=1,MCHAM2.NOMCHE(/2)
  382.             IF(MCHAM2.NOMCHE(I_COMP).EQ. MOCOMP) THEN
  383.               IF (MCHAM2.TYPCHE(I_COMP).NE.'REAL*8') THEN
  384.                 MOTERR(1:16)  = MCHAM2.TYPCHE(I_COMP)
  385.                 MOTERR(17:20) = MOCOMP
  386.                 MOTERR(21:36) = MCHEL2.TITCHE
  387.                 CALL ERREUR(552)
  388.                 RETURN
  389.               ENDIF
  390.               MELVA2= MCHAM2.IELVAL(I_COMP)
  391.               NPT2  = MELVA2.VELCHE(/1)
  392.               NEL2  = MELVA2.VELCHE(/2)
  393.               GOTO 471
  394.             ENDIF
  395.   461     CONTINUE
  396.           MOTERR = MOCOMP
  397.           CALL ERREUR(677)
  398.           RETURN
  399.   471     CONTINUE
  400.  
  401.           DO I_EL1=1,NB_EL1
  402.             q_totE = XZERO
  403.             q_dejE = XZERO
  404. C           1ere boucle pour connaitre le cumul dans l'element
  405.             DO I_NO1=1,NB_NO1
  406.               x_prop = MELVA1.VELCHE(MIN(I_NO1,NPT1),MIN(I_EL1,NEL1))
  407.               x_qtot = MELVA2.VELCHE(MIN(I_NO1,NPT2),MIN(I_EL1,NEL2))
  408.               q_totE = q_totE + x_qtot
  409.               q_dejE = q_dejE + x_qtot*x_prop
  410.             ENDDO
  411.  
  412.             DO I_NO1=1,NB_NO1
  413.               NU_INC = IPT1.NUM(I_NO1,I_EL1)
  414.               S_ELEM = ICPR4(NU_INC,I_CONS)
  415.               NU_EL  = S_ELEM.NB_EL + 1
  416.               S_ELEM.NB_EL        = NU_EL
  417.               S_ELEM.MELPRO(NU_EL)= MELVA1
  418.               S_ELEM.LISNEL(NU_EL)= I_EL1
  419.               S_ELEM.LISNNO(NU_EL)= I_NO1
  420.  
  421.               x_prop = MELVA1.VELCHE(MIN(I_NO1,NPT1),MIN(I_EL1,NEL1))
  422.               x_qtot = MELVA2.VELCHE(MIN(I_NO1,NPT2),MIN(I_EL1,NEL2))
  423.               x_deja = x_prop*x_qtot
  424.               S_ELEM.QTOTNO(NU_EL)= x_qtot
  425.               S_ELEM.QDEJNO(NU_EL)= x_deja
  426.               S_ELEM.QRESNO(NU_EL)= x_qtot - x_deja
  427.  
  428.               S_ELEM.QTOTEL(NU_EL)= q_totE
  429.               S_ELEM.QDEJEL(NU_EL)= q_dejE
  430.               S_ELEM.QRESEL(NU_EL)= q_totE - q_dejE
  431.             ENDDO
  432.           ENDDO
  433.         ENDDO
  434.       ENDDO
  435.  
  436. * +--------------------------------------------------------------------+
  437. * |    Remplissage du S_LX a partir des 'LX' du MCHPOI                 |
  438. * +--------------------------------------------------------------------+
  439.       DO 30 I_LX=1,N
  440.         NUM_LX = MELEME.NUM(1,I_LX)
  441.         NU_INC = SCPR.ICPR1(NUM_LX)
  442.         IF(NU_INC .EQ. 0)GOTO 30
  443.         I_CONS = SCPR.ICPR2(NUM_LX)
  444.  
  445.         S_LX   = S_CONS.LIS_LX(I_CONS)
  446.         N_LX   = S_LX.NB_LX + 1
  447.         S_ELEM = ICPR4(NU_INC,I_CONS)
  448.         qmax=XZERO
  449.         DO I_EL=1,S_ELEM.NB_EL
  450.           qmax = qmax + S_ELEM.QTOTNO(I_EL)
  451.         ENDDO
  452.         S_LX.NB_LX       = N_LX
  453.         S_LX.LISELE(N_LX)= S_ELEM
  454.         S_LX.LISNOL(N_LX)= NUM_LX
  455.         S_LX.LISNOI(N_LX)= NU_INC
  456.         S_LX.Q_REA(N_LX) = MPOVAL.VPOCHA(I_LX,1)
  457.         S_LX.Q_MAX(N_LX) = qmax
  458.  30   CONTINUE
  459.  
  460. * +--------------------------------------------------------------------+
  461. * |    Boucle pour repartir la chaleur sur les noeuds des elements     |
  462. * +--------------------------------------------------------------------+
  463.       SEGINI,ITRAV
  464.       ISEGSU(**)=ITRAV
  465.       DO 60 I_CONS=1,NB_CON
  466.         S_LX=S_CONS.LIS_LX(I_CONS)
  467.         DO 70 I_LX=1,S_LX.NB_LX
  468.           S_ELEM = S_LX.LISELE(I_LX)
  469.           NUM_LX = S_LX.LISNOL(I_LX)
  470.           NU_INC = S_LX.LISNOI(I_LX)
  471.           N_EL   = S_ELEM.NB_EL
  472.           qrea   = S_LX.Q_REA(I_LX)
  473.           qmax   = S_LX.Q_MAX(I_LX)
  474.           qprec  = MAX(ABS(qmax)*XZPREL, XPETIT/XZPREC)
  475.  
  476.           IF(qrea .GT. XZERO)THEN
  477. C           AU CHAUFFAGE
  478. C           Premier Trie selon QDEJEL pour connaitre les elements qui ont commence la transformation
  479.             DO I_EL=1,N_EL
  480.               ITRAV.IPERM1(I_EL)=I_EL
  481.               ITRAV.TT2(I_EL)   =-S_ELEM.QDEJEL(I_EL)
  482.             ENDDO
  483.             CALL TRIFLO(ITRAV.TT2,ITRAV.TT1,
  484.      &                  ITRAV.IPERM1,ITRAV.IT1,N_EL)
  485.  
  486. C           On compte combien ont commence
  487.             NB_COM=0
  488.             DO I_EL=1,N_EL
  489.               iperm = ITRAV.IPERM1(I_EL)
  490.               q_totE= S_ELEM.QTOTEL(iperm)
  491.               q_dejE= S_ELEM.QDEJEL(iperm)
  492.               IF (q_dejE .GT. (XZPREL*q_totE)) THEN
  493.                 NB_COM=NB_COM + 1
  494.                 ITRAV.TT2(NB_COM)=S_ELEM.QRESNO(iperm)
  495.               ENDIF
  496.             ENDDO
  497.             NB_NUL=N_EL-NB_COM
  498.  
  499.             IF(NB_COM .GT. 1)THEN
  500. C             Deuxieme Trie : Les elements qui ont commences sont tries selon QRESNO
  501.               CALL TRIFLO(ITRAV.TT2,ITRAV.TT1,
  502.      &                    ITRAV.IPERM1,ITRAV.IT1,NB_COM)
  503.             ENDIF
  504.  
  505.             IF(NB_NUL .GT. 1)THEN
  506. C             Troisieme Trie : Les elements qui n'ont pas commences sont tries selon QTOTNO
  507.               DO I_EL=1,NB_NUL
  508.                 iperm          = ITRAV.IPERM1(NB_COM+I_EL)
  509.                 ITRAV.TT2(I_EL)= S_ELEM.QTOTNO(iperm)
  510.               ENDDO
  511.               CALL TRIFLO(ITRAV.TT2,ITRAV.TT1,
  512.      &                    ITRAV.IPERM1(NB_COM+1),ITRAV.IT1,NB_NUL)
  513.             ENDIF
  514.  
  515. C           Il ne reste qu'a repartir le plus equitablement possible aux noeuds
  516. C           PRINT * ,'Hausse ELEMENT',N_EL,NB_COM,NB_NUL,qrea,S_ELEM
  517.             DO 80 I_EL=1,N_EL
  518.               iperm =ITRAV.IPERM1(I_EL)
  519.               melvap=S_ELEM.MELPRO(iperm)
  520.               q_tot =S_ELEM.QTOTNO(iperm)
  521.               q_dej =S_ELEM.QDEJNO(iperm)
  522.               q_res =S_ELEM.QRESNO(iperm)
  523.               IP    =S_ELEM.LISNNO(iperm)
  524.               IE    =S_ELEM.LISNEL(iperm)
  525.  
  526.               IF(NB_COM .GE. 1)THEN
  527.                 q_part = qrea / NB_COM
  528.               ELSE
  529.                 q_part = qrea / NB_NUL
  530.               ENDIF
  531.  
  532.               IF(q_tot.EQ.XZERO .OR. (q_res-q_part) .LE. XZERO)THEN
  533. C               Ce noeud termine sa transformation vers le haut
  534.                 MELVAP.VELCHE(IP,IE)=XUN
  535.                 qrea  = qrea - q_res
  536. C               PRINT *,'Hausse-a',qrea,q_tot,q_dej,q_res,q_part
  537.                 q_res = XZERO
  538.                 q_dej = q_tot
  539.  
  540.               ELSE
  541. C               Ce noeud consomme exactement q_part
  542.                 pp    = MAX(MIN((q_dej+q_part)/q_tot,XUN),XZERO)
  543.                 MELVAP.VELCHE(IP,IE)=pp
  544.                 qrea  = qrea  - q_part
  545. C               PRINT *,'Hausse-b',qrea,q_tot,q_dej,q_res,q_part
  546.                 q_res = q_res - q_part
  547.                 q_dej = q_tot - q_res
  548.               ENDIF
  549.  
  550.               NB_COM = NB_COM - 1
  551.               IF(NB_COM.LT.0)THEN
  552.                 NB_NUL = NB_NUL - 1
  553.               ENDIF
  554.  
  555.               S_ELEM.QRESNO(iperm)=q_res
  556.               S_ELEM.QDEJNO(iperm)=q_dej
  557.  80         CONTINUE
  558.             IF(ABS(qrea).GT.qprec) CALL soucis(1)
  559.  
  560.           ELSEIF(qrea .LT. XZERO)THEN
  561. C           AU REFROIDISSEMENT
  562. C           Premier Trie selon QRESEL pour connaitre les elements qui ont commence la transformation inverse
  563.             DO I_EL=1,N_EL
  564.               ITRAV.IPERM1(I_EL)= I_EL
  565.               ITRAV.TT2(I_EL)   =-S_ELEM.QRESEL(I_EL)
  566.             ENDDO
  567.             CALL TRIFLO(ITRAV.TT2,ITRAV.TT1,
  568.      &                  ITRAV.IPERM1,ITRAV.IT1,N_EL)
  569.  
  570. C           On compte combien d'elements ont commence la transformation inverse
  571.             NB_COM=0
  572.             DO I_EL=1,N_EL
  573.               iperm = ITRAV.IPERM1(I_EL)
  574.               q_totE= S_ELEM.QTOTEL(iperm)
  575.               q_resE= S_ELEM.QRESEL(iperm)
  576.               IF (q_resE .GT. (XZPREL*q_totE)) THEN
  577.                 NB_COM=NB_COM + 1
  578.                 ITRAV.TT2(NB_COM)=S_ELEM.QDEJNO(iperm)
  579.               ENDIF
  580.             ENDDO
  581.             NB_NUL=N_EL-NB_COM
  582.  
  583.             IF(NB_COM .GT. 1)THEN
  584. C             Deuxieme Trie : Les elements qui ont commences sont tries selon QDEJNO
  585.               CALL TRIFLO(ITRAV.TT2,ITRAV.TT1,
  586.      &                    ITRAV.IPERM1,ITRAV.IT1,NB_COM)
  587.             ENDIF
  588.  
  589.             IF(NB_NUL .GT. 1)THEN
  590. C             Troisieme Trie : Les elements qui n'ont pas commences sont tries selon QTOTNO
  591.               DO I_EL=1,NB_NUL
  592.                 iperm          = ITRAV.IPERM1(NB_COM+I_EL)
  593.                 ITRAV.TT2(I_EL)= S_ELEM.QTOTNO(iperm)
  594.               ENDDO
  595.               CALL TRIFLO(ITRAV.TT2,ITRAV.TT1,
  596.      &                    ITRAV.IPERM1(NB_COM+1),ITRAV.IT1,NB_NUL)
  597.             ENDIF
  598.  
  599. C           Il ne reste qu'a repartir le plus equitablement possible aux noeuds
  600. C           PRINT * ,'Baisse ELEMENT',N_EL,NB_COM,NB_NUL,qrea
  601.             DO 90 I_EL=1,N_EL
  602.               iperm =ITRAV.IPERM1(I_EL)
  603.               melvap=S_ELEM.MELPRO(iperm)
  604.               q_tot =S_ELEM.QTOTNO(iperm)
  605.               q_dej =S_ELEM.QDEJNO(iperm)
  606.               q_res =S_ELEM.QRESNO(iperm)
  607.               IP    =S_ELEM.LISNNO(iperm)
  608.               IE    =S_ELEM.LISNEL(iperm)
  609.  
  610.               IF(NB_COM .GE. 1)THEN
  611.                 q_part = qrea / NB_COM
  612.               ELSE
  613.                 q_part = qrea / NB_NUL
  614.               ENDIF
  615.  
  616.               IF(q_tot.EQ.XZERO .OR. (q_dej+q_part) .LE. XZERO)THEN
  617. C               Ce noeud termine sa transformation vers le bas
  618.                 MELVAP.VELCHE(IP,IE)=XZERO
  619.                 qrea  = qrea + q_dej
  620. C               PRINT *,'Baisse-a',qrea,q_tot,q_dej,q_res,q_part
  621.                 q_res = q_tot
  622.                 q_dej = XZERO
  623.  
  624.               ELSE
  625. C               Ce noeud consomme exactement q_part
  626.                 pp    = MAX(MIN((q_dej+q_part)/q_tot,XUN),XZERO)
  627.                 MELVAP.VELCHE(IP,IE)=pp
  628.                 qrea  = qrea  - q_part
  629. C               PRINT *,'Baisse-b',qrea,q_tot,q_dej,q_res,q_part
  630.                 q_res = q_res - q_part
  631.                 q_dej = q_tot - q_res
  632.               ENDIF
  633.  
  634.               NB_COM = NB_COM - 1
  635.               IF(NB_COM.LT.0)THEN
  636.                 NB_NUL = NB_NUL - 1
  637.               ENDIF
  638.  
  639.               S_ELEM.QRESNO(iperm)=q_res
  640.               S_ELEM.QDEJNO(iperm)=q_dej
  641.  90         CONTINUE
  642.             IF(ABS(qrea).GT.qprec) CALL soucis(1)
  643.  
  644.           ELSE
  645. C           Cas (qrea .EQ. XZERO) : survient si chaleur latente nulle
  646. C           On fait donc une bascule des proportions de phase
  647. C           Si pp>0 ==> pp=1
  648. C           Si pp<1 ==> pp=0
  649.             DO 100 I_EL=1,N_EL
  650.               MELVAP=S_ELEM.MELPRO(I_EL)
  651.               IP    =S_ELEM.LISNNO(I_EL)
  652.               IE    =S_ELEM.LISNEL(I_EL)
  653.               pp    =MELVAP.VELCHE(IP,IE)
  654.               IF    (pp .LT. XUN  )THEN
  655.                 pp = XUN
  656.               ELSEIF(pp .GT. XZERO)THEN
  657.                 pp = XZERO
  658.               ELSE
  659.                 CALL ERREUR(5)
  660.               ENDIF
  661.               MELVAP.VELCHE(IP,IE)=pp
  662.  100        CONTINUE
  663.           ENDIF
  664.  70     CONTINUE
  665.  60   CONTINUE
  666.  
  667. C     On supprime les SEGMENTS de travail
  668.       ISEGSU(**)=ISEGSU
  669.       DO ii=1,ISEGSU(/1)
  670.         ISEG=ISEGSU(ii)
  671.         SEGSUP,ISEG
  672.       ENDDO
  673.  
  674.       CALL ACTOBJ('MCHAML',MCHELP,1)
  675.       CALL ECROBJ('MCHAML',MCHELP)
  676.       END
  677.  
  678.  
  679.  
  680.  
  681.  
  682.  

© Cast3M 2003 - Tous droits réservés.
Mentions légales