Télécharger la plaquette Cast3M


fr en

Télécharger ycls.eso

Retour à la liste

Numérotation des lignes : 

ycls

  1. C YCLS      SOURCE    CB215821  24/04/12    21:17:31     11897          
  2.       SUBROUTINE YCLS(NOMPR,MTABX,MTABZ,IHV,MZIN,NOMI,MATRIK,MCHPO1)
  3.       IMPLICIT INTEGER(I-N)
  4.       IMPLICIT REAL*8 (A-H,O-Z)
  5. C***********************************************************************
  6. C
  7. C  CET OPERATEUR DISCRETISE Convection Diffusion + source
  8. C
  9. C  SYNTAXE :
  10. C  ---------
  11. C
  12. C
  13. C  COEFFICIENTS :
  14. C  --------------
  15. C
  16. C
  17. C  ALF  (SCAL DOMA)  DIFFUSIVITE THERMIQUE
  18. C       (SCAL ELEM)
  19. C       (SCAL NOEU)
  20. C
  21. C  INCONNUES :
  22. C  -----------
  23. C
  24. C   TN                 CHAMP DE TEMPERATURE
  25. C   UN                 CHAMP DE VITESSE TRANSPORTANT
  26. C
  27. C   NOMPR Nom de l'opérateur appelant (NS,TSCA,KONV ou LAPN)
  28. C         suivant l'opérateur les coefs sont différents
  29. C   Schéma en temps ou l'opérateur intervient :
  30. C    SEMI            CNG                         TVISQ
  31. C    teta - schéma   Crank Nicolson généralisé   Tenseur visqueux
  32. C   Ces schéma nécessitent que DFDT soit en CENTREE
  33. C   Le Schéma BDF2 ne fait intervenir que DFDT
  34. C
  35. C E/  MTABX Table de l'opérateur -> coefficients
  36. C E/  MTABZ Table DOMAINE fourni par la lecture du modèle (.dgibi)
  37. C E/  IHV=0 SCALAIRE
  38. C     IHV=1 VITESSE
  39. C E/  MZIN  pointeur du CHPOINT de l'inconnue
  40. C E/  NOMI  nom de l'inconnue
  41. C************************************************************************
  42.  
  43.  
  44. -INC PPARAM
  45. -INC CCOPTIO
  46. -INC CCREEL
  47. -INC CCGEOME                                                            
  48. -INC SMCOORD
  49. -INC SMLENTI
  50.       POINTEUR MLENT4.MLENTI
  51. -INC SMLMOTS
  52. -INC SMMODEL
  53. -INC SMELEME
  54.       POINTEUR IGEOM.MELEME,MELEMS.MELEME,MELEMC.MELEME,MELEMP.MELEME
  55.       POINTEUR MELEM2.MELEME
  56. -INC SMCHPOI
  57.       POINTEUR MZIN.MCHPOI
  58.       POINTEUR MTETA1.MPOVAL,MTETA2.MPOVAL,MTETA3.MPOVAL,MTETA4.MPOVAL
  59. -INC SMTABLE
  60.       POINTEUR MTABZ.MTABLE
  61. -INC SIZFFB
  62.       POINTEUR IZF1.IZFFM,IZH2.IZHR,IZW.IZFFM,IZWH.IZHR
  63.       SEGMENT SAJT
  64.       REAL*8 AJT(IDIM,IDIM,NPG),RF1(NP,MP,IDIM),SM1(NP,IDIM)
  65.       REAL*8 TN1(NP,IDIM),TN2(NP,IDIM)
  66.       REAL*8 WT(MP,NPG), WS(MP,NPG)
  67.       REAL*8 WTC(IDIM*NPG,MP),WSC(IDIM*NPG,MP),HRD(IDIM*NPG,NP)
  68.       REAL*8 FND(NPG,NP) ,FNE(NPG,NP,IDIM),HRE(NPG,NP,IDIM)
  69.       REAL*8 GMU(NPG,IDIM),USI(IDIM)
  70.       ENDSEGMENT
  71. -INC SMRIGID
  72. C-INC SMMATRIK
  73. C*******************************************************************
  74. C
  75.       SEGMENT MATRIK
  76.       REAL*8  COEMTK(NMATRI)
  77.       INTEGER JRIGEL(NRIGE,NMATRI)
  78.       INTEGER KSYM,KMINC,KMINCP,KMINCD,KIZM
  79.       INTEGER KISPGT,KISPGP,KISPGD
  80.       INTEGER KNTTT,KNTTP,KNTTD
  81.       INTEGER KIDMAT(NKID)
  82.       INTEGER KKMMT(NKMT)
  83.       ENDSEGMENT
  84.  
  85.       SEGMENT JMATRI
  86.       CHARACTER*(LOCOMP) LISPRI(NBMF),LJSDUA(NBMF)
  87.       INTEGER LIZAFM(NBSOUS,NBMF)
  88.       INTEGER KSPGP,KSPGD
  89.       ENDSEGMENT
  90.       POINTEUR JMATRS.JMATRI,JMATR1.JMATRI,JMATR2.JMATRI,JMATR3.JMATRI
  91.  
  92. C Stokage matrices elementaires non assemblees (valeurs)
  93.       SEGMENT IZAFM
  94.       REAL*8  AM(NBEL,NP,MP)
  95.       ENDSEGMENT
  96.       POINTEUR IPM1.IZAFM,IPM2.IZAFM,IPM3.IZAFM,IPM4.IZAFM
  97.  
  98. C*******************************************************************
  99.       POINTEUR IPM.IZAFM
  100. -INC SMCHAML
  101.       POINTEUR MCHELG.MCHELM,MCHAMG.MCHAML,MELVAG.MELVAL
  102.  
  103.       LOGICAL TDFDT,TCONV,TSOUR,TSOURB,TMATRI
  104.       LOGICAL ICAL,XPG,XRIG
  105.       LOGICAL XDIAG,XTV,XTG,XBDF,XCONS
  106.       CHARACTER*8 CHHH,TYPE,NOM,NOM0,CHAI,SCHT,NOMPER,MPRE
  107.       CHARACTER*4 NMACO
  108.       CHARACTER*(*) NOMI,NOMPR
  109.       CHARACTER*(LOCOMP) NOMP(3),NOMD(3)
  110.       DIMENSION ICOEF(4),SU(3),XPOI(3)
  111. C*****************************************************************************
  112. CTCLSF
  113.       ICAL=.FALSE.
  114.  
  115.       NOMPER=NOMPR
  116. c       write(6,*)' DEBUT YCLS appele par ',NOMPER
  117.  
  118. C- Table de l'opérateur (pointeur MTABX)
  119.  
  120. C- Récupération de la table RV (pointeur MTAB1)
  121.       CALL LEKTAB(MTABX,'EQEX',MTAB1)
  122. C- Récupération de la table des options de l'opérateur (pointeur KOPTI)
  123.       CALL LEKTAB(MTABX,'KOPT',KOPTI)
  124. C- Récupération de la table DOMAINE de la zone (pointeur MTABZ) ABANDONNE
  125. C- Pour le parallélisme on récupère la table DOMAINE du modèle partitionné donné
  126. C- dans Gibiane (pointeur MTABZ)
  127. C      CALL LEKTAB(MTABX,'DOMZ',MTABZ)   ABANDONNE
  128. C- Récupération de la table INCO (pointeur KINC)
  129.       CALL LEKTAB(MTAB1,'INCO',KINC)
  130.       IF(IERR.NE.0) RETURN
  131.  
  132. C*****************************************************************************
  133. C Traitement des options
  134. C KIMPL = 0 -> EXPL     1 -> IMPL     2 -> CN
  135. C                                   AIMPL=0 explicite
  136. C IKOMP = 0 Formulation non conservative des termes de convection
  137. C IKOMP = 1 Formulation faible conservative des termes de convection
  138. C           i.e. integre par partie on fait apparaitre l'integrale de surface du
  139. C           flux convectif ce qui permet d'imposer le cas echeant le flux total
  140. C           (convectif + diffusif)
  141. C IKOMP = 2 Correction de Temam pour la stabilite L2
  142. C KMU = 0 QDM : mu = cste sur l'element
  143. C KMU = 1 QDM : Formulation en Tau (contrainte visqueuse)
  144. C KMU = 2 QDM : Formulation en Grad de Mu
  145. C KKALE=1  Indice indiquant que l'on est en ALE
  146. C sinon KKALE=0
  147. C KFORM = 0 -> SI       1 -> EF       2 -> VF    3 -> EFMC
  148. C IDCEN  : entier indiquant le type de décentrement souhaité
  149. C IDCEN 1-> CENTREE  2-> SUPGDC  3-> SUPG  4-> TVISQUEU  5-> CNG
  150. C E/  CMD    : Coefficient multiplicateur du décentrement
  151. C Si IDCEN=4 ou =5 CMD=DT
  152.  
  153.       CALL ACME(KOPTI,'KFORM',KFORM)
  154.       CALL ACME(KOPTI,'KIMPL',KIMPL)
  155.       CALL ACMF(KOPTI,'AIMPL',AIMPL)
  156.       CALL ACME(KOPTI,'IKOMP',IKOMP)
  157.       CALL ACME(KOPTI,'KMACO',KMACO)
  158.       KKMACO=KMACO
  159.       CALL ACME(KOPTI,'IDCEN',IDCEN)
  160.       CALL ACMF(KOPTI,'CMD',CMD)
  161.       CALL ACME(KOPTI,'ALE',KKALE)
  162.       CALL ACME(KOPTI,'KMU',KMU)
  163.       IF(KMU.EQ.2.AND.IHV.EQ.1)THEN
  164.       IHM=1
  165.       ELSE
  166.       IHM=0
  167.       ENDIF
  168.  
  169. C Restrictions des options
  170.       IF(NOMPER.NE.'NS')KKALE=0
  171.       IF(NOMPER.EQ.'LAPN')IDCEN=1
  172.  
  173.       IF(KIMPL.EQ.0)AIMPL=0.D0
  174.       IF(KIMPL.EQ.1)AIMPL=1.D0
  175.  
  176. c       write(6,*)' YCLS KMACO=',KMACO,'IDCEN=',IDCEN
  177.        IF(KMACO.EQ.1.AND.IDCEN.NE.5)THEN
  178.          CALL ACMM(KOPTI,'NMACO',NMACO)
  179.          TYPE=' '
  180.          CALL ACMO(MTAB1,NMACO,TYPE,MATRIK)
  181.         IF(TYPE.EQ.'MATRIK')THEN
  182.           KKMACO=2
  183.         ELSE
  184.           KKMACO=1
  185.         ENDIF
  186.        ENDIF
  187.       IF(IDCEN.EQ.5)KKMACO=0
  188. c      write(6,*)' YCLS KMACO=',kmaco,' KKMACO=',KKMACO,'IKOMP=',ikomp
  189.  
  190. C*****************************************************************************
  191.  
  192.  
  193.  
  194.       KPOIND=0
  195.       NBMF=1
  196.       XRIG=.FALSE.
  197.       TMATRI=.TRUE.
  198. C      C2LAPN=1.D0 Laplacien   C2LAPN=0.D0 pas de laplacien
  199.       C2LAPN=1.D0
  200.       TCONV=.TRUE.
  201.       TSOUR=.FALSE.
  202.       TSOURB=.FALSE.
  203.       XPG=.FALSE.
  204.       IF(IDCEN.GT.1)XPG=.TRUE.
  205.       IF(IDCEN.EQ.4)AIMPL=0.D0
  206. C      C1CNG cas particulier de Crank Nicolson généralisé
  207.        C1CNG=1.D0
  208.       IF(IDCEN.EQ.5)THEN
  209.        AIMPL=1.D0
  210.        C1CNG=0.5D0
  211.       ENDIF
  212.  
  213.        AIEX=AIMPL
  214.       IF(AIMPL.EQ.0.D0.AND.KKMACO.NE.0)THEN
  215.        TMATRI=.TRUE.
  216.        AIEX=1.D0
  217.       ENDIF
  218.  
  219.       IF(AIMPL.EQ.0.D0.AND.KKMACO.EQ.0)THEN
  220.        TMATRI=.FALSE.
  221.       ENDIF
  222.  
  223.       IAXI=0
  224.       IF(IFOMOD.EQ.0)IAXI=2
  225.       DEUPI=1.D0
  226.       IF(IAXI.NE.0)DEUPI=2.D0*XPI
  227.  
  228.  
  229.         IF(XRIG)THEN
  230.          IF(IHV.EQ.0)THEN
  231.           NOMP(1)='T   '
  232.           NOMD(1)='Q   '
  233.          ELSEIF(IHV.EQ.1)THEN
  234.           NOMP(1)='UX  '
  235.           NOMP(2)='UY  '
  236.           NOMP(3)='UZ  '
  237.           NOMD(1)='FX  '
  238.           NOMD(2)='FY  '
  239.           NOMD(3)='FZ  '
  240.          ENDIF
  241.         ELSE
  242.          IF(IHV.EQ.0)THEN
  243.           NOMP(1)=NOMI
  244.           NOMD(1)=NOMI
  245.          ELSEIF(IHV.EQ.1)THEN
  246.           WRITE(NOMP(1),FMT='(I1)')1
  247.           WRITE(NOMD(1),FMT='(I1)')1
  248.           WRITE(NOMP(2),FMT='(I1)')2
  249.           WRITE(NOMD(2),FMT='(I1)')2
  250.           WRITE(NOMP(3),FMT='(I1)')3
  251.           WRITE(NOMD(3),FMT='(I1)')3
  252.           NOMP(1)=NOMP(1)(1:1)//NOMI(1:LOCOMP-1)
  253.           NOMD(1)=NOMD(1)(1:1)//NOMI(1:LOCOMP-1)
  254.           NOMP(2)=NOMP(2)(1:1)//NOMI(1:LOCOMP-1)
  255.           NOMD(2)=NOMD(2)(1:1)//NOMI(1:LOCOMP-1)
  256.           NOMP(3)=NOMP(3)(1:1)//NOMI(1:LOCOMP-1)
  257.           NOMD(3)=NOMD(3)(1:1)//NOMI(1:LOCOMP-1)
  258.  
  259.          ENDIF
  260.         ENDIF
  261. c       write(6,*)' NOMP=',nomp
  262. c       write(6,*)' NOMD=',nomd
  263.  
  264.       CALL LEKTAB(MTABZ,'SOMMET',MELEMS)
  265.       CALL LEKTAB(MTABZ,'MAILLAGE',MELEME)
  266.       MELEM2=MELEME
  267.  
  268. C  On cree un second membre non vide
  269.             SEGACT MELEMS
  270.             N=MELEMS.NUM(/2)
  271.             IF(IHV.EQ.0)NC=1
  272.             IF(IHV.EQ.1)NC=IDIM
  273.             SEGDES MELEMS
  274.             NSOUPO=1
  275.             NAT=2
  276.             SEGINI MCHPO1,MSOUP1,MPOVA1
  277.             MCHPO1.JATTRI(1)=2
  278.             MCHPO1.IFOPOI=IFOUR
  279.             MCHPO1.MTYPOI='    '
  280.             MCHPO1.MOCHDE='                                            '
  281.             MCHPO1.IPCHP(1)=MSOUP1
  282.             MSOUP1.IGEOC=MELEMS
  283.             MSOUP1.IPOVAL=MPOVA1
  284.              IF(IHV.EQ.0)THEN
  285.               MSOUP1.NOCOMP(1)=NOMD(1)
  286.              ELSE
  287.               DO 91 N=1,NC
  288.               MSOUP1.NOCOMP(N)=NOMD(N)
  289.  91           CONTINUE
  290.              ENDIF
  291.             SEGDES MSOUP1,MCHPO1
  292.  
  293. C*************************************************************************
  294. C Lecture des Inco(s) aux temps precedents si Transitoire
  295. c       write(6,*)'Lecture des INCO '
  296.  
  297.         CALL LRCHT(MZIN,MTETA1,TYPE,IGEOM)
  298.         CALL KRIPAD(IGEOM,MLENT1)
  299.         CALL VERPAD(MLENT1,MELEMS,IRET)
  300.          IF(IRET.NE.0)THEN
  301. C%    Indice %m1:8 : L'objet %m9:16 n'a pas le bon support géométrique
  302.           MOTERR(1: 8) = 'TETA1'
  303.           MOTERR(9:16) = 'CHPOINT '
  304.           WRITE(IOIMP,*)'Opérateur : ',NOMPER
  305.           CALL ERREUR(788)
  306.           RETURN
  307.          ENDIF
  308.  
  309.         NC=MTETA1.VPOCHA(/2)
  310.  
  311.         IF((IHV.EQ.0.AND.NC.NE.1).OR.
  312.      &     (IHV.EQ.1.AND.NC.NE.IDIM))THEN
  313. C% Indice %m1:8 : L'objet %m9:16 n'a pas le bon nombre de composantes
  314.          MOTERR( 1: 8) = 'Inconnue'
  315.          MOTERR( 9:16) = 'CHPOINT'
  316.          WRITE(IOIMP,*)'Operateur : ',NOMPER
  317.          CALL ERREUR(784)
  318.          RETURN
  319.         ENDIF
  320.  
  321.         SEGDES IGEOM,MELEMS
  322.  
  323. C*****************************************************************************
  324. C Lecture du coefficient
  325. c       write(6,*)'Lecture des COEF '
  326.       NMUVARI=0
  327.       MCHELG =0
  328.  
  329.       CALL ACME(MTABX,'IARG',IARG)
  330. c      write(6,*)' YCLS IARG=',IARG,' KKALE=',KKALE
  331.  
  332. C--------- CAS NS ou TSCA Formulation non conservative ou conservative -------
  333.        IF((NOMPER.EQ.'NS'.OR.NOMPER.EQ.'TSCA').AND.
  334.      & (IKOMP.EQ.0.OR.IKOMP.EQ.1))THEN
  335.  
  336. c_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._
  337.        IF( (IHV.EQ.0.AND.(IARG.EQ.3.OR.IARG.EQ.4)).OR.
  338.      c     (IHV.EQ.1.AND.(IARG.EQ.3.OR.IARG.EQ.4.OR.IARG.EQ.6))
  339.      c                                            )THEN
  340.        C2LAPN=1.D0
  341.        TCONV=.TRUE.
  342.        TSOUR=.FALSE.
  343.        TSOURB=.FALSE.
  344.        IF(IARG.EQ.4)TSOUR=.TRUE.
  345.        IF(IHV.EQ.1.AND.IARG.EQ.6)THEN
  346.        TSOUR=.TRUE.
  347.        TSOURB=.TRUE.
  348.        ENDIF
  349. C Ro
  350.        NUCOEF=1
  351.        CALL LEKMOF(MTABZ,NUCOEF,MTABX,KINC,0,MCHEL1,KPOIND,0,MUVARI)
  352.        IF (IERR.NE.0) RETURN
  353. C Un
  354.         NUCOEF=NUCOEF+1
  355.         CALL LEKMOF(MTABZ,NUCOEF,MTABX,KINC,1,MCHEL3,KPOIND,0,MUVARI)
  356.         IF (IERR.NE.0) RETURN
  357. C Mu
  358.        NUCOEF=NUCOEF+1
  359.        CALL LEKMOF(MTABZ,NUCOEF,MTABX,KINC,0,MCHEL2,KPOIND,IHM,MUVARI)
  360.        IF(MUVARI.NE.0)THEN
  361.        MCHELG=MUVARI
  362.        NMUVARI=NMUVARI+1
  363.        ENDIF
  364.        IF (IERR.NE.0) RETURN
  365.  
  366.        ELSE
  367.       WRITE(6,*)' Nombre d''arguments incorrect : ',IARG
  368.       IF(IHV.EQ.0)THEN
  369.       WRITE(6,*)' On attend 3 ou 4'
  370.       ENDIF
  371.       IF(IHV.EQ.1)THEN
  372.       WRITE(6,*)' On attend 3 ou 6'
  373.       ENDIF
  374. C           Indice %m1:8 : nombre d'arguments incorrect
  375.             MOTERR(1:8) = NOMPER
  376.             CALL ERREUR(804)
  377.          RETURN
  378.        ENDIF
  379. c_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._
  380.  
  381.  
  382.        IF(TSOUR)THEN
  383. C Source -> MCHEL4
  384. C Gb
  385.         NUCOEF=NUCOEF+1
  386. c      write(6,*)' Coef 3 ',NUCOEF,' MTABZ=',MTABZ
  387.         IVC=IHV
  388.         CALL LEKMOF(MTABZ,NUCOEF,MTABX,KINC,IVC,MCHEL4,KPOIND,0,MUVARI)
  389.         IF (IERR.NE.0) RETURN
  390.  
  391.         IF(TSOURB)THEN
  392. C Source Boussinesq pour la QDM : Gb(T0 - T) -> MCHEL4
  393. C Tn
  394. c      write(6,*)' Coef 4 BOU',NUCOEF
  395.          NUCOEF=NUCOEF+1
  396.          CALL LEKMOF(MTABZ,NUCOEF,MTABX,KINC,0,MCHEL5,KPOIND,0,MUVARI)
  397.          IF (IERR.NE.0) RETURN
  398. C To
  399.          NUCOEF=NUCOEF+1
  400.          CALL LEKMOF(MTABZ,NUCOEF,MTABX,KINC,0,MCHEL6,KPOIND,0,MUVARI)
  401.          IF (IERR.NE.0) RETURN
  402. C -> MCHEL4 = Gb*(Tn - T0)
  403.          CALL MELBOU(MTABZ,MCHEL4,MCHEL4,MCHEL5,MCHEL6)
  404.         ENDIF
  405.        ENDIF
  406.  
  407. C       write(6,*)'FIN CAS NS ou TSCA Formulation non conservative'
  408. C----- FIN CAS NS ou TSCA Formulation non conservative ou non conservative ---
  409.  
  410. C--------- CAS NS ou TSCA Formulation Temam ----------------------------------
  411.        ELSEIF((                  NOMPER.EQ.'TSCA').AND.IKOMP.EQ.2)THEN
  412.  
  413.        TCONV=.TRUE.
  414.        TSOUR=.TRUE.
  415.        TSOURB=.FALSE.
  416. C On convient d'avoir 3 coefs tout le temps (mu , u , s)
  417.  
  418. C Nu
  419.        NUCOEF=1
  420. c      write(6,*)' Coef 1 ',NUCOEF,' MTABZ=',MTABZ
  421.        CALL LEKMOF(MTABZ,NUCOEF,MTABX,KINC,0,MCHEL2,KPOIND,IHM,MUVARI)
  422.        IF(MUVARI.NE.0)THEN
  423.        MCHELG=MUVARI
  424.        NMUVARI=NMUVARI+1
  425.        ENDIF
  426.        IF (IERR.NE.0) RETURN
  427.  
  428.        SEGACT MELEME
  429.        L1=72
  430.        N1=MAX(1,LISOUS(/1))
  431.        N2=1
  432.        N3=6
  433.        SEGINI MCHEL1
  434.        SEGACT MCHEL2
  435.        DO 75 L=1,N1
  436.        MCHAM2=MCHEL2.ICHAML(L)
  437.        SEGACT MCHAM2
  438.        MELVA2=MCHAM2.IELVAL(1)
  439.        SEGACT MELVA2
  440.        N1PTEL=MELVA2.VELCHE(/1)
  441.        N1EL=1
  442.        N2PTEL=0
  443.        N2EL=0
  444.        SEGINI MCHAM1,MELVA1
  445.        MCHEL1.ICHAML(L)=MCHAM1
  446.        MCHAM1.IELVAL(1)=MELVA1
  447.        DO 76 LG=1,N1PTEL
  448.        MELVA1.VELCHE(LG,1)=1.D0
  449.  76    CONTINUE
  450.  75    CONTINUE
  451.  
  452. C Un
  453.         NUCOEF=NUCOEF+1
  454. c      write(6,*)' Coef 2 ',NUCOEF
  455.         CALL LEKMOF(MTABZ,NUCOEF,MTABX,KINC,1,MCHEL3,KPOIND,0,MUVARI)
  456.         IF (IERR.NE.0) RETURN
  457.         CALL LEKDIV(MTABZ,NUCOEF,MTABX,KINC,MCHEL5,KPOIND)
  458.         IF (IERR.NE.0) RETURN
  459.  
  460.        IF(TSOUR)THEN
  461. C Source -> MCHEL4
  462. C Gb
  463.         NUCOEF=NUCOEF+1
  464. c      write(6,*)' Coef 3 ',NUCOEF,' MTABZ=',MTABZ
  465.         IVC=IHV
  466.         CALL LEKMOF(MTABZ,NUCOEF,MTABX,KINC,IVC,MCHEL4,KPOIND,0,MUVARI)
  467.         IF (IERR.NE.0) RETURN
  468.        ENDIF
  469.  
  470. C----- FIN CAS NS ou TSCA Formulation Temam ----------------------------------
  471.  
  472. C--------- CAS LAPN ----------------------------------------------------------
  473.        ELSEIF(NOMPER.EQ.'LAPN')THEN
  474.  
  475.        TSOUR=.FALSE.
  476.        TSOURB=.FALSE.
  477.        TCONV =.FALSE.
  478. C Alfa
  479.        NUCOEF=1
  480. c      write(6,*)' Coef 1 ',NUCOEF
  481.        CALL LEKMOF(MTABZ,NUCOEF,MTABX,KINC,0,MCHEL2,KPOIND,IHM,MUVARI)
  482.        IF(MUVARI.NE.0)THEN
  483.        MCHELG=MUVARI
  484.        NMUVARI=NMUVARI+1
  485.        ENDIF
  486.        IF (IERR.NE.0) RETURN
  487.  
  488.        SEGACT MELEME
  489.        L1=72
  490.        N1=MAX(1,LISOUS(/1))
  491.        N2=1
  492.        N3=6
  493.        SEGINI MCHEL1,MCHEL3
  494.        SEGACT MCHEL2
  495.        DO 74 L=1,N1
  496.        MCHAM2=MCHEL2.ICHAML(L)
  497.        SEGACT MCHAM2
  498.        MELVA2=MCHAM2.IELVAL(1)
  499.        SEGACT MELVA2
  500.        N1PTEL=MELVA2.VELCHE(/1)
  501.        N1EL=1
  502.        N2PTEL=0
  503.        N2EL=0
  504.        SEGINI MCHAM1,MELVA1
  505.        MCHEL1.ICHAML(L)=MCHAM1
  506.        MCHAM1.IELVAL(1)=MELVA1
  507.        DO 73 LG=1,N1PTEL
  508.        MELVA1.VELCHE(LG,1)=0.D0
  509.  73    CONTINUE
  510.        N1PTEL=N1PTEL*IDIM
  511.        SEGINI MCHAM3,MELVA3
  512.        MCHEL3.ICHAML(L)=MCHAM3
  513.        MCHAM3.IELVAL(1)=MELVA3
  514.        DO 731 LG=1,N1PTEL
  515.        MELVA3.VELCHE(LG,1)=0.D0
  516.  731   CONTINUE
  517.  74    CONTINUE
  518. C----- FIN CAS LAPN ----------------------------------------------------------
  519.  
  520. C--------- CAS KONV ----------------------------------------------------------
  521.        ELSEIF(NOMPER.EQ.'KONV')THEN
  522.  
  523.        IF(IARG.NE.3.AND.IARG.NE.4)THEN
  524.       WRITE(6,*)' Nombre d''arguments incorrect : ',IARG
  525.       WRITE(6,*)' On attend 3 ou 4'
  526. C           Indice %m1:8 : nombre d'arguments incorrect
  527.             MOTERR(1:8) = NOMPER
  528.             CALL ERREUR(804)
  529.          RETURN
  530.       ENDIF
  531.  
  532.        C2LAPN=0.D0
  533.        TCONV=.TRUE.
  534.        TSOUR=.FALSE.
  535.        TSOURB=.FALSE.
  536. C Ro
  537.        NUCOEF=1
  538. c      write(6,*)' Coef 1 ',NUCOEF
  539.        CALL LEKMOF(MTABZ,NUCOEF,MTABX,KINC,0,MCHEL1,KPOIND,0,MUVARI)
  540.        IF (IERR.NE.0) RETURN
  541.  
  542. C Un
  543.         NUCOEF=NUCOEF+1
  544. c      write(6,*)' Coef 2 ',NUCOEF
  545.         CALL LEKMOF(MTABZ,NUCOEF,MTABX,KINC,1,MCHEL3,KPOIND,0,MUVARI)
  546.         IF (IERR.NE.0) RETURN
  547. C Mu
  548.        NUCOEF=NUCOEF+1
  549. c      write(6,*)' Coef 3 ',NUCOEF
  550.        CALL LEKMOF(MTABZ,NUCOEF,MTABX,KINC,0,MCHEL2,KPOIND,0,MUVARI)
  551.        IF (IERR.NE.0) RETURN
  552.  
  553. C Dt
  554.        IF(IDCEN.EQ.4.OR.IDCEN.EQ.5)THEN
  555.        NUCOEF=NUCOEF+1
  556. c      write(6,*)' Coef 4 ',NUCOEF
  557.        CALL LEKMOF(MTABZ,NUCOEF,MTABX,KINC,0,MCHEL4,KPOIND,0,MUVARI)
  558.        IF (IERR.NE.0) RETURN
  559.       SEGACT MCHEL4
  560.       MCHAM4=MCHEL4.ICHAML(1)
  561.       SEGACT MCHAM4
  562.       MELVA4=MCHAM4.IELVAL(1)
  563.       SEGACT MELVA4
  564.       N1PTEL=MELVA4.VELCHE(/1)
  565.       N1EL=MELVA4.VELCHE(/2)
  566.         IF(N1EL.EQ.1)THEN
  567.         DT= MELVA4.VELCHE(1,1)
  568.         SEGSUP MELVA4,MCHAM4,MCHEL4
  569.         ELSE
  570.         call erreur(991)
  571.         RETURN
  572.         ENDIF
  573.  
  574.        CMD=DT
  575.        ENDIF
  576.  
  577. C----- FIN CAS KONV ----------------------------------------------------------
  578.  
  579.        ENDIF
  580. c      write(6,*)' Fin lect COEF '
  581.         IF (IERR.NE.0) RETURN
  582. C*************************************************************************
  583. C******* CALCUL **********************************************************
  584. C*************************************************************************
  585.  
  586.        SEGACT MTETA1
  587.  
  588.       CALL KRIPAD(MELEMS,MLENTI)
  589.       SEGDES MELEMS
  590.  
  591.       SEGACT MELEME,MELEM2
  592.       NBSOUS=LISOUS(/1)
  593.       IF(NBSOUS.EQ.0)NBSOUS=1
  594.  
  595. C MATRIK . MATRIK . MATRIK . MATRIK . MATRIK . MATRIK  MATRIK . MATRIK .
  596.       IF(KKMACO.NE.2)THEN
  597. C ... Création / Ajustement des matrices élémentaires
  598.       IF(TMATRI)THEN
  599. C Cas RIGIDITE
  600.        IF(XRIG)THEN
  601.  
  602.          NRIGE = 8
  603.          NRI   = 0
  604.          NRIGEL=NRI+NBSOUS
  605.          SEGINI MRIGID
  606.          MRIGID.MTYMAT = '        '
  607.          MRIGID.IFORIG = IFOUR
  608.  
  609. c     write(6,*)' NRIGE,NRIGEL,MRIGID=',NRIGE,NRIGEL,MRIGID
  610.  
  611.        ELSE
  612. C Cas MATRIK
  613.  
  614.          NRIGE=7
  615.          NKID =9
  616.          NKMT =7
  617.          NMATRI=0
  618.          SEGINI MATRIK
  619. c        write(6,*)' YCLS creation MATRIK =',matrik
  620.  
  621.          NRIGE=7
  622.          NKID =9
  623.          NKMT =7
  624.          NMATR0=JRIGEL(/2)
  625.  
  626.          NBME=1
  627.          IF(IHV.EQ.1)NBME=IDIM
  628.          NMATRI=NBME
  629.          SEGADJ MATRIK
  630.  
  631.          DO 41 M=1,NBME
  632.          JRIGEL(1,M)=MELEME
  633.          JRIGEL(2,M)=MELEM2
  634.          JRIGEL(7,M)=0
  635.          IF(TCONV)JRIGEL(7,M)=2
  636.          SEGINI JMATRI
  637.          JRIGEL(4,M)=JMATRI
  638.          KSPGP=MELEMS
  639.          KSPGD=MELEMS
  640.          LISPRI(1)=NOMP(M)//'    '
  641.          LJSDUA(1)=NOMD(M)//'    '
  642.  41      CONTINUE
  643.        ENDIF
  644.       ENDIF
  645. C ... Fin Création / Ajustement des matrices élémentaires
  646.  
  647.  
  648.       ELSE
  649. C Activation RIGIDITE
  650.        IF(XRIG)THEN
  651.          SEGACT MRIGID
  652.        ELSE
  653. C Activation MATRIK
  654.          SEGACT MATRIK
  655.        ENDIF
  656.  
  657. C ... Activation des matrices élémentaires
  658.       IF(TMATRI)THEN
  659. C Cas RIGIDITE
  660.        IF(XRIG)THEN
  661.          NRIGE=8
  662.          NRI=IRIGEL(/2)
  663.          NRIGEL=NBSOUS+NRI
  664.          SEGAct MRIGID
  665. c     write(6,*)' NRIGE,NRIGEL,MRIGID=',NRIGE,NRIGEL,MRIGID
  666.  
  667.        ELSE
  668. C Cas MATRIK Pleine
  669.          NMATRI=JRIGEL(/2)
  670.          NBME=NMATRI
  671.  
  672.          DO 42  M=1,NBME
  673.          JMATRI=JRIGEL(4,M)
  674.          SEGACT JMATRI*MOD
  675.          LISPRI(1)=NOMP(M)//'    '
  676.          LJSDUA(1)=NOMD(M)//'    '
  677.  42      CONTINUE
  678.  
  679.        ENDIF
  680.       ENDIF
  681. C ... Fin Activation des matrices élémentaires
  682.  
  683.       ENDIF
  684. C FIN MATRIK . FIN MATRIK . FIN MATRIK . FIN MATRIK . FIN MATRIK . FIN M
  685. C ......................................................................
  686.  
  687.       SEGACT MCHEL1
  688.       SEGACT MCHEL2
  689.       IF(MCHELG.NE.0)THEN
  690.         IF(NMUVARI.GT.1)THEN
  691.         WRITE(6,*)' Trop de coefficients pour la diffusion ?!?'
  692. C% Données incompatibles
  693.         CALL ERREUR(22)
  694.         RETURN
  695.         ENDIF
  696.       SEGACT MCHELG
  697.       ENDIF
  698.       SEGACT MCHEL3
  699.       IF(TSOUR)THEN
  700.       SEGACT MCHEL4
  701.       ENDIF
  702.       IF(IKOMP.EQ.2)THEN
  703.       SEGACT MCHEL5
  704.       ENDIF
  705.  
  706.        IF(MAX(1,MELEM2.LISOUS(/1)).NE.MAX(1,LISOUS(/1)))THEN
  707.         WRITE(6,*)' Geometries incompatibles dans ',nomper
  708. C% Données incompatibles
  709.         CALL ERREUR(22)
  710.         RETURN
  711.         ENDIF
  712.  
  713.       NKD=0
  714.       DO 101 L=1,MAX(1,LISOUS(/1))
  715.       IPT1=MELEME
  716.       IPT2=MELEM2
  717.       IF(LISOUS(/1).NE.0)IPT1=LISOUS(L)
  718.       IF(MELEM2.LISOUS(/1).NE.0)THEN
  719.       IPT2=MELEM2.LISOUS(L)
  720.       NKD=0
  721.       ENDIF
  722.       SEGACT IPT1,IPT2
  723.  
  724.       NOM0 = NOMS(IPT1.ITYPEL)//'    '
  725.       CALL KALPBG(NOM0,'FONFORM ',IZFFM)
  726.       IF(IZFFM.EQ.0)RETURN
  727.       SEGACT IZFFM*MOD
  728.       IZHR=KZHR(1)
  729.       SEGACT IZHR*MOD
  730.       IZF1 = KTP(1)
  731.       IZH2 = KZHR(2)
  732.       IZW  = IZFFM
  733.       IZWH = IZHR
  734.  
  735.       NES=GR(/1)
  736.       NPG=GR(/3)
  737.  
  738.       NP = IPT1.NUM(/1)
  739.       MP = IPT2.NUM(/1)
  740. C?    MP = IZW.FN(/1) ceci doit etre identique
  741.  
  742.       NBEL=IPT1.NUM(/2)
  743.  
  744.       SEGINI SAJT
  745.  
  746. c.......................................................................
  747. C MATRIK . MATRIK . MATRIK . MATRIK . MATRIK . MATRIK  MATRIK . MATRIK .
  748. c     write(6,*)' KKMACO=',KKMACO,'TMATRI=',TMATRI
  749.       IF(KKMACO.NE.2)THEN
  750.        IF(TMATRI)THEN
  751. C Cas RIGIDITE
  752.         IF(XRIG)THEN
  753.         IRIGEL(1,NRI+L)=MELEME
  754.         COERIG(L)=1.D0
  755.  
  756.         IRIGEL(7,NRI+L)=0
  757.         IF(TCONV)IRIGEL(7,NRI+L)=2
  758.  
  759.         NBME=1
  760.         IF(IHV.EQ.1)NBME=IDIM
  761.         NLIGRP=NP
  762.         NLIGRD=MP
  763.         SEGINI DESCR
  764.         IRIGEL(3,NRI+L)=DESCR
  765.          IF(NBME.EQ.1)THEN
  766.            DO 102 I=1,NLIGRP
  767.             LISINC(I)=NOMP(1)//'    '
  768.             NOELEP(I)=I
  769.  102       CONTINUE
  770.            DO 103 I=1,NLIGRD
  771.             LISDUA(I)=NOMD(1)//'    '
  772.             NOELED(I)=I
  773.  103       CONTINUE
  774.          ELSE
  775.          ENDIF
  776.         SEGDES DESCR
  777.  
  778.         NELRIG=NBEL
  779.         SEGINI xMATRI
  780.         IRIGEL(4,NRI+L)=xMATRI
  781.         xmatri.symre=irigel(7,nri+l)
  782. c       write(6,*)'NELRIG,IMATRI=',NELRIG,IMATRI
  783.  
  784. *        DO 104 K=1,NELRIG
  785. *        SEGINI XMATRI
  786. c       write(6,*)'NLIGRD,NLIGRP,XMATRI=',NLIGRD,NLIGRP,XMATRI
  787. *        IMATTT(K)=XMATRI
  788. * 104    CONTINUE
  789.  
  790.         ELSE
  791. C Cas MATRIK
  792.         SEGINI IPM1
  793.         mtail=(IPM1.AM(/1))*(ipm1.am(/2))*(ipm1.am(/3))
  794.         JMATR1=JRIGEL(4,NMATR0+1)
  795.         JMATR1.LIZAFM(L,1)=IPM1
  796.         IPM2=IPM1
  797.         IPM3=IPM1
  798.          IF(NBME.GE.2)THEN
  799.           JMATR2=JRIGEL(4,NMATR0+2)
  800.           IF(IAXI.NE.0)THEN
  801.            SEGINI IPM2
  802.         mtail=(IPM2.AM(/1))*(ipm2.am(/2))*(ipm2.am(/3))
  803.            JMATR2.LIZAFM(L,1)=IPM2
  804.            ICAL=.TRUE.
  805.           ELSE
  806.            IPM2=IPM1
  807.            JMATR2.LIZAFM(L,1)=IPM2
  808.            ICAL=.FALSE.
  809.           ENDIF
  810.          ENDIF
  811.          IF(NBME.GE.3)THEN
  812.           JMATR3=JRIGEL(4,NMATR0+3)
  813.            IPM3=IPM1
  814.            JMATR3.LIZAFM(L,1)=IPM3
  815.            ICAL=.FALSE.
  816.          ENDIF
  817.        ENDIF
  818.       ENDIF
  819. c       write(6,*)' YCLS appelé par ',NOMPER,' Taille matrice=',mtail
  820.  
  821.       ELSE
  822.  
  823. C Cas RIGIDITE
  824.        IF(XRIG)THEN
  825.  
  826.        ELSE
  827. C Cas MATRIK
  828.         JMATR1=JRIGEL(4,1)
  829.         SEGACT JMATR1
  830.         IPM1=JMATR1.LIZAFM(L,1)
  831.         IPM2=IPM1
  832.         IPM3=IPM1
  833.         SEGACT IPM1
  834.          IF(NBME.GE.2)THEN
  835.           JMATR2=JRIGEL(4,2)
  836.           SEGACT JMATR2
  837.           IPM2=JMATR2.LIZAFM(L,1)
  838.           SEGACT IPM2
  839.          ENDIF
  840.          IF(NBME.GE.3)THEN
  841.           JMATR3=JRIGEL(4,3)
  842.           IPM3=JMATR3.LIZAFM(L,1)
  843.           SEGACT IPM3
  844.          ENDIF
  845.        ENDIF
  846.       ENDIF
  847. C FIN MATRIK . FIN MATRIK . FIN MATRIK . FIN MATRIK . FIN MATRIK . FIN M
  848. C ......................................................................
  849.  
  850.  
  851. C----Ro = 1.
  852.       IK1=1
  853.       MCHAM1=MCHEL1.ICHAML(L)
  854.       SEGACT MCHAM1
  855.  
  856.       MELVA1=MCHAM1.IELVAL(1)
  857.       SEGACT MELVA1
  858.       N1PTEL=MELVA1.VELCHE(/1)
  859.       N1EL=MELVA1.VELCHE(/2)
  860.         IF(N1EL.EQ.1)THEN
  861.         IK1=1
  862.         ELSEIF(N1EL.EQ.NBEL)THEN
  863.         IK1=0
  864.         ENDIF
  865.  
  866. C----Nu ou Lambda
  867.       IK2=1
  868.       MCHAM2=MCHEL2.ICHAML(L)
  869.       SEGACT MCHAM2
  870.  
  871.       MELVA2=MCHAM2.IELVAL(1)
  872.       SEGACT MELVA2
  873.       N1PTEL=MELVA2.VELCHE(/1)
  874.       N1EL=MELVA2.VELCHE(/2)
  875.         IF(N1EL.EQ.1)THEN
  876.         IK2=1
  877.         ELSEIF(N1EL.EQ.NBEL)THEN
  878.         IK2=0
  879.         ENDIF
  880.  
  881.       IF(MCHELG.NE.0)THEN
  882.       MCHAMG=MCHELG.ICHAML(L)
  883.       SEGACT MCHAMG
  884.  
  885.       MELVAG=MCHAMG.IELVAL(1)
  886.       SEGACT MELVAG
  887.       ENDIF
  888.  
  889. C----U
  890.       IK3=1
  891.       MCHAM3=MCHEL3.ICHAML(L)
  892.       SEGACT MCHAM3
  893.  
  894.       MELVA3=MCHAM3.IELVAL(1)
  895.       SEGACT MELVA3
  896.       N1PTEL=MELVA3.VELCHE(/1)
  897.       N1EL=MELVA3.VELCHE(/2)
  898.         IF(N1EL.EQ.1)THEN
  899.         IK3=1
  900.         ELSEIF(N1EL.EQ.NBEL)THEN
  901.         IK3=0
  902.         ENDIF
  903.  
  904. C----sources
  905.       IK4=1
  906.        IF(TSOUR)THEN
  907.       MCHAM4=MCHEL4.ICHAML(L)
  908.       SEGACT MCHAM4
  909.  
  910.       MELVA4=MCHAM4.IELVAL(1)
  911.       SEGACT MELVA4
  912.       N1PTEL=MELVA4.VELCHE(/1)
  913.       N1EL=MELVA4.VELCHE(/2)
  914.         IF(N1EL.EQ.1)THEN
  915.         IK4=1
  916.         ELSEIF(N1EL.EQ.NBEL)THEN
  917.         IK4=0
  918.         ENDIF
  919.        ENDIF
  920.  
  921. C----Div (rocp u )
  922.       IK5=1
  923.        IF(IKOMP.EQ.2)THEN
  924.       MCHAM5=MCHEL5.ICHAML(L)
  925.       SEGACT MCHAM5
  926.  
  927.       MELVA5=MCHAM5.IELVAL(1)
  928.       SEGACT MELVA5
  929.       N1PTEL=MELVA5.VELCHE(/1)
  930.       N1EL=MELVA5.VELCHE(/2)
  931.         IF(N1EL.EQ.1)THEN
  932.         IK5=1
  933.         ELSEIF(N1EL.EQ.NBEL)THEN
  934.         IK5=0
  935.         ENDIF
  936.        ENDIF
  937.  
  938. C      write(6,*)' AVANT 108 NC=',NC,' NBEL=',NBEL,MP,NP,NC
  939. C===============================================
  940.         AI1=AIMPL-1.D0
  941.         IF(AIMPL.EQ.0.D0)THEN
  942.           AI2=-1.D0
  943.         ELSE
  944.           AI2=AI1/AIMPL
  945.         ENDIF
  946.  
  947.         DO 108 KE=1,NBEL
  948.  
  949.         NK1=KE + IK1*(1 - KE)
  950.         NK2=KE + IK2*(1 - KE)
  951.         NK3=KE + IK3*(1 - KE)
  952.         NK4=KE + IK4*(1 - KE)
  953.         NK5=KE + IK5*(1 - KE)
  954.  
  955.         DO     I=1,NP
  956.         J=IPT1.NUM(I,KE)
  957.         DO     N=1,IDIM
  958.         XYZ(N,I)=XCOOR((J-1)*(IDIM+1)+N)
  959.         ENDDO
  960.         ENDDO
  961.  
  962.         DO     I=1,NP
  963.          I1=MLENT1.LECT(IPT1.NUM(I,KE))
  964.         DO     N=1,NC
  965.          TN1(I,N)=MTETA1.VPOCHA(I1,N)
  966.         ENDDO
  967.         ENDDO
  968.  
  969.         CALL CALJBR(FN,GR,PG,XYZ,HR,PGSQ,RPG,NES,
  970.      *  IDIM,NP,NPG,IAXI,AIRE,AJ,SGN)
  971.  
  972.         IF (IDCEN.EQ.0.OR.IDCEN.EQ.1) THEN
  973.         CALL RSETD(WT,FN,(NP*NPG))
  974.         CALL RSETD(WS,FN,(NP*NPG))
  975.         ELSE
  976.         CALL CALDEC(WT,WS,XYZ,GR,HR,FN,NES,IDIM,NP,NPG,AJT,
  977.      &  IDCEN,CMD,MELVA1.VELCHE(1,NK1),MELVA2.VELCHE(1,NK2),
  978.      &  MELVA3.VELCHE(1,NK3),TN1,NC,IKOMP,XREF,AIRE,KE)
  979.         ENDIF
  980.  
  981.         CALL INITD(SM1,(MP*IDIM),0.D0)
  982.  
  983. C:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
  984.         IF(KKMACO.EQ.2)THEN
  985.  
  986. C ...... Chargement Rigidite ou Matrik
  987. C Cas RIGIDITE
  988.           IF(XRIG)THEN
  989. *           XMATRI=IMATTT(KE)
  990.            DO  I=1,MP
  991.            DO  J=1,NP
  992.            RE(I,J,ke)=RF1(J,I,1)
  993.            ENDDO
  994.            ENDDO
  995. *           SEGDES XMATRI
  996.           ELSE
  997. C Cas MATRIK
  998.            DO  N=1,NBME
  999.            JMATR1=JRIGEL(4,N)
  1000.            IPM4=JMATR1.LIZAFM(L,1)
  1001.            DO  I=1,MP
  1002.            DO  J=1,NP
  1003.            RF1(J,I,N)=IPM4.AM(KE,J,I)
  1004.            ENDDO
  1005.            ENDDO
  1006.            ENDDO
  1007.           ENDIF
  1008. C...... Multiplication Matrice Vecteur
  1009.  
  1010.         DO 710 I=1,MP
  1011.         DO 711 J=1,NP
  1012.  
  1013.         DO 716 N=1,NC
  1014.           SM1(I,N)   = SM1(I,N) + AI2*RF1(J,I,N)*TN1(J,N)
  1015.  716    CONTINUE
  1016.  
  1017.  711    CONTINUE
  1018.  710    CONTINUE
  1019.  
  1020.         ELSEIF(KKMACO.NE.2)THEN
  1021. C.................................................................
  1022. C...... Convection  Laplacien ....  préparation forme conservative
  1023. C...... Convection  Laplacien ....             et non conservative
  1024. C préparation
  1025.           AKOMP=1.D0
  1026.           IF(IKOMP.EQ.1)AKOMP=-1.D0
  1027.           NIPG=NPG*IDIM
  1028.           I1=0
  1029.         DO  LG=1,NPG
  1030.           PDR=PGSQ(LG)*DEUPI*RPG(LG)
  1031.           C1=C1CNG*MELVA1.VELCHE(LG,NK1)*PDR*AKOMP
  1032.           C2=C2LAPN*MELVA2.VELCHE(LG,NK2)*PDR
  1033.          DO  N=1,IDIM
  1034.           I1=I1+1
  1035.           C3=MELVA3.VELCHE((N-1)*NPG+LG,NK3)*C1
  1036.          DO  I=1,MP
  1037.           WTC(I1,I)=WT(I,LG)*C3 + IZWH.HR(N,I,LG)*C2
  1038.          ENDDO
  1039.          ENDDO
  1040.         ENDDO
  1041.  
  1042.           I1=0
  1043.         DO     LG=1,NPG
  1044.          DO     N=1,IDIM
  1045.           I1=I1+1
  1046.          DO     J=1,NP
  1047.           HRD(I1,J)=HR(N,J,LG)
  1048.          ENDDO
  1049.          ENDDO
  1050.         ENDDO
  1051.  
  1052.        IF(IKOMP.EQ.0.OR.IKOMP.EQ.2)THEN
  1053. C...... Convection .... forme non conservative
  1054.         DO     I=1,MP
  1055.         DO     J=1,NP
  1056.           U3=0.D0
  1057.         DO 402 I1=1,NIPG
  1058.           U3=U3+(WTC(I1,I) * HRD(I1,J))
  1059.  402      CONTINUE
  1060.  
  1061.         DO 403 N=1,NC
  1062.           SM1(I,N)   = SM1(I,N) + AI1*U3*TN1(J,N)
  1063.  403    CONTINUE
  1064.           RF1(J,I,1) = AIEX*U3
  1065.  
  1066.         ENDDO
  1067.         ENDDO
  1068.        ENDIF
  1069.  
  1070.        IF(IKOMP.EQ.1)THEN
  1071. C...... Convection .... forme conservative
  1072.         DO     I=1,MP
  1073.         DO     J=1,NP
  1074.           U3=0.D0
  1075.         DO 412 I1=1,NIPG
  1076.           U3=U3+(WTC(I1,I) * HRD(I1,J))
  1077.  412      CONTINUE
  1078.  
  1079.         DO 413 N=1,NC
  1080.           SM1(J,N)   = SM1(J,N) + AI1*U3*TN1(I,N)
  1081.  413    CONTINUE
  1082.           RF1(I,J,1) = AIEX*U3
  1083.  
  1084.         ENDDO
  1085.         ENDDO
  1086.        ENDIF
  1087.  
  1088. C...... Fin Convection ....
  1089.  
  1090. C.......... Cas KMU=1 et FTAU pour QDM (Formulation en To)
  1091.       IF(KMU.EQ.1)THEN
  1092. c       write(6,*)' Cas FTAU'
  1093.  
  1094.         DO      LG=1,NPG
  1095.           PDR=PGSQ(LG)*DEUPI*RPG(LG)
  1096.           C2=C2LAPN*MELVA2.VELCHE(LG,NK2)*PDR
  1097.          DO      N=1,IDIM
  1098.          DO      J=1,NP
  1099.           FNE(LG,J,N)=HR(N,J,LG)*C2
  1100.           HRE(LG,J,N)=HR(N,J,LG)
  1101.          ENDDO
  1102.          ENDDO
  1103.         ENDDO
  1104.  
  1105.       IF(IDIM.EQ.2)THEN
  1106.         DO 4192 I=1,NP
  1107.           US1=0.D0
  1108.           US2=0.D0
  1109.         DO 4182 J=1,NP
  1110.           U11=0.D0
  1111.           U12=0.D0
  1112.           U21=0.D0
  1113.           U22=0.D0
  1114.         DO 4172 LG=1,NPG
  1115.           U11=U11+HRE(LG,J,1)*FNE(LG,I,1)
  1116.           U12=U12+HRE(LG,J,1)*FNE(LG,I,2)
  1117.           U21=U21+HRE(LG,J,2)*FNE(LG,I,1)
  1118.           U22=U22+HRE(LG,J,2)*FNE(LG,I,2)
  1119.  4172   CONTINUE
  1120.           US1 = US1 - U11*TN1(J,1) - U12*TN1(J,2)
  1121.           US2 = US2 - U21*TN1(J,1) - U22*TN1(J,2)
  1122.  4182   CONTINUE
  1123.           SM1(I,1) = SM1(I,1) + US1
  1124.           SM1(I,2) = SM1(I,2) + US2
  1125.  4192   CONTINUE
  1126.       ELSE
  1127.         DO 5193 I=1,NP
  1128.           US1=0.D0
  1129.           US2=0.D0
  1130.           US3=0.D0
  1131.         DO 5183 J=1,NP
  1132.           U11=0.D0
  1133.           U12=0.D0
  1134.           U13=0.D0
  1135.           U21=0.D0
  1136.           U22=0.D0
  1137.           U23=0.D0
  1138.           U31=0.D0
  1139.           U32=0.D0
  1140.           U33=0.D0
  1141.         DO 5173 LG=1,NPG
  1142.           U11=U11+HRE(LG,J,1)*FNE(LG,I,1)
  1143.           U12=U12+HRE(LG,J,1)*FNE(LG,I,2)
  1144.           U13=U13+HRE(LG,J,1)*FNE(LG,I,3)
  1145.  
  1146.           U21=U21+HRE(LG,J,2)*FNE(LG,I,1)
  1147.           U22=U22+HRE(LG,J,2)*FNE(LG,I,2)
  1148.           U23=U23+HRE(LG,J,2)*FNE(LG,I,3)
  1149.  
  1150.           U31=U31+HRE(LG,J,3)*FNE(LG,I,1)
  1151.           U32=U32+HRE(LG,J,3)*FNE(LG,I,2)
  1152.           U33=U33+HRE(LG,J,3)*FNE(LG,I,3)
  1153.  5173   CONTINUE
  1154.           US1 = US1 - U11*TN1(J,1) - U12*TN1(J,2) - U13*TN1(J,3)
  1155.           US2 = US2 - U21*TN1(J,1) - U22*TN1(J,2) - U23*TN1(J,3)
  1156.           US3 = US3 - U31*TN1(J,1) - U32*TN1(J,2) - U33*TN1(J,3)
  1157.  5183   CONTINUE
  1158.           SM1(I,1) = SM1(I,1) + US1
  1159.           SM1(I,2) = SM1(I,2) + US2
  1160.           SM1(I,3) = SM1(I,3) + US3
  1161.  5193   CONTINUE
  1162.       ENDIF
  1163.  
  1164.  
  1165.       ENDIF
  1166. C...... FIN Cas KMU=1 et FTAU pour QDM (Formulation en To)
  1167.  
  1168. C.......... Cas KMU=2 et MUVARI pour QDM (Formulation en Grad de Mu)
  1169.       IF(KMU.EQ.2.AND.MCHELG.NE.0)THEN
  1170. c       write(6,*)' Cas MUVARI'
  1171.  
  1172.           CALL INITD(GMU,(NPG*IDIM),0.D0)
  1173.         DO     LG=1,NPG
  1174.         DO     I=1,NP
  1175.         DO     N=1,NC
  1176.           GMU(LG,N)=GMU(LG,N)+HR(N,I,LG)*MELVAG.VELCHE(I,NK2)
  1177.         ENDDO
  1178.         ENDDO
  1179.         ENDDO
  1180.  
  1181.  
  1182.         DO     LG=1,NPG
  1183.           PDR=PGSQ(LG)*DEUPI*RPG(LG)
  1184.          DO     N=1,IDIM
  1185.          DO     J=1,NP
  1186.           FNE(LG,J,N)=FN(J,LG)*GMU(LG,N)
  1187.           HRE(LG,J,N)=HR(N,J,LG)*PDR
  1188.          ENDDO
  1189.          ENDDO
  1190.          ENDDO
  1191.  
  1192.       IF(IDIM.EQ.2)THEN
  1193.         DO 419 I=1,NP
  1194.           US1=0.D0
  1195.           US2=0.D0
  1196.         DO 418 J=1,NP
  1197.           U11=0.D0
  1198.           U12=0.D0
  1199.           U21=0.D0
  1200.           U22=0.D0
  1201.         DO 417 LG=1,NPG
  1202.           U11=U11+HRE(LG,J,1)*FNE(LG,I,1)
  1203.           U12=U12+HRE(LG,J,1)*FNE(LG,I,2)
  1204.           U21=U21+HRE(LG,J,2)*FNE(LG,I,1)
  1205.           U22=U22+HRE(LG,J,2)*FNE(LG,I,2)
  1206.  417    CONTINUE
  1207.           US1 = US1 - U11*TN1(J,1) - U12*TN1(J,2)
  1208.           US2 = US2 - U21*TN1(J,1) - U22*TN1(J,2)
  1209.  418    CONTINUE
  1210.           SM1(I,1) = SM1(I,1) + US1
  1211.           SM1(I,2) = SM1(I,2) + US2
  1212.  419    CONTINUE
  1213.       ELSE
  1214.         DO 4193 I=1,NP
  1215.           US1=0.D0
  1216.           US2=0.D0
  1217.           US3=0.D0
  1218.         DO 4183 J=1,NP
  1219.           U11=0.D0
  1220.           U12=0.D0
  1221.           U13=0.D0
  1222.           U21=0.D0
  1223.           U22=0.D0
  1224.           U23=0.D0
  1225.           U31=0.D0
  1226.           U32=0.D0
  1227.           U33=0.D0
  1228.         DO 4173 LG=1,NPG
  1229. c         U11=U11+FN(I,LG)*HRE(LG,J,1)*GMU(LG,1)
  1230. c         U12=U12+FN(I,LG)*HRE(LG,J,1)*GMU(LG,2)
  1231. c         U21=U21+FN(I,LG)*HRE(LG,J,2)*GMU(LG,1)
  1232. c         U22=U22+FN(I,LG)*HRE(LG,J,2)*GMU(LG,2)
  1233.  
  1234.           U11=U11+HRE(LG,J,1)*FNE(LG,I,1)
  1235.           U12=U12+HRE(LG,J,1)*FNE(LG,I,2)
  1236.           U13=U13+HRE(LG,J,1)*FNE(LG,I,3)
  1237.  
  1238.           U21=U21+HRE(LG,J,2)*FNE(LG,I,1)
  1239.           U22=U22+HRE(LG,J,2)*FNE(LG,I,2)
  1240.           U23=U23+HRE(LG,J,2)*FNE(LG,I,3)
  1241.  
  1242.           U31=U31+HRE(LG,J,3)*FNE(LG,I,1)
  1243.           U32=U32+HRE(LG,J,3)*FNE(LG,I,2)
  1244.           U33=U33+HRE(LG,J,3)*FNE(LG,I,3)
  1245.  4173   CONTINUE
  1246.           US1 = US1 - U11*TN1(J,1) - U12*TN1(J,2) - U13*TN1(J,3)
  1247.           US2 = US2 - U21*TN1(J,1) - U22*TN1(J,2) - U23*TN1(J,3)
  1248.           US3 = US3 - U31*TN1(J,1) - U32*TN1(J,2) - U33*TN1(J,3)
  1249.  4183   CONTINUE
  1250.           SM1(I,1) = SM1(I,1) + US1
  1251.           SM1(I,2) = SM1(I,2) + US2
  1252.           SM1(I,3) = SM1(I,3) + US3
  1253.  4193   CONTINUE
  1254.       ENDIF
  1255.  
  1256.  
  1257.       ENDIF
  1258. C...... FIN Cas KMU=2 et MUVARI pour QDM (Formulation en Grad de Mu)
  1259.  
  1260. C.......... Cas CNG uniquement (complément)
  1261.         IF(IDCEN.EQ.5)THEN
  1262. C préparation
  1263.           I1=0
  1264.         DO     LG=1,NPG
  1265.           PDR=PGSQ(LG)*DEUPI*RPG(LG)
  1266.           C1=C1CNG*MELVA1.VELCHE(LG,NK1)*PDR
  1267.          DO     N=1,IDIM
  1268.           I1=I1+1
  1269.           C3=MELVA3.VELCHE((N-1)*NPG+LG,NK3)*C1
  1270.          DO     I=1,MP
  1271.           WSC(I1,I)=WS(I,LG)*C3
  1272.          ENDDO
  1273.          ENDDO
  1274.         ENDDO
  1275.  
  1276.         DO 510 I=1,MP
  1277.         DO 511 J=1,NP
  1278.           U3=0.D0
  1279.         DO 512 I1=1,NIPG
  1280.           U3=U3+WSC(I1,I) * HRD(I1,J)
  1281.  512      CONTINUE
  1282.  
  1283.         DO 516 N=1,NC
  1284.           SM1(I,N)   = SM1(I,N) - U3*TN1(J,N)
  1285.  516    CONTINUE
  1286.  
  1287.  511    CONTINUE
  1288.  510    CONTINUE
  1289.         ENDIF
  1290. C..........Fin Cas CNG uniquement (complément)
  1291.  
  1292. C..............................................................
  1293. C..............................................................
  1294.         IF(IKOMP.EQ.2)THEN
  1295. C...... Convection  Laplacien .... Formulation conservative
  1296. C préparation
  1297.         DO     LG=1,NPG
  1298.           PDR=PGSQ(LG)*DEUPI*RPG(LG)
  1299.           C5=C1CNG*MELVA5.VELCHE(LG,NK5)*PDR
  1300.          DO     I=1,MP
  1301.           WTC(LG,I)=WT(I,LG)*C5
  1302.          ENDDO
  1303.         ENDDO
  1304.  
  1305.         DO     LG=1,NPG
  1306.          DO     J=1,NP
  1307.           FND(LG,J)=FN(J,LG)
  1308.          ENDDO
  1309.         ENDDO
  1310.  
  1311.         DO 440 I=1,MP
  1312.         DO 441 J=1,NP
  1313.           U3=0.D0
  1314.         DO 442 LG=1,NPG
  1315.           U3=U3+WTC(LG,I) * FND(LG,J)
  1316.  442      CONTINUE
  1317.  
  1318.         DO 446 N=1,NC
  1319.           SM1(I,N)   = SM1(I,N) + AI1*U3*TN1(J,N)
  1320.  446    CONTINUE
  1321.           RF1(J,I,1)=RF1(J,I,1)+(AIEX*U3)
  1322.  
  1323.  441    CONTINUE
  1324.  440    CONTINUE
  1325.  
  1326. C.......... Cas CNG uniquement (complément)
  1327.         IF(IDCEN.EQ.5)THEN
  1328. C préparation
  1329.           I1=0
  1330.         DO     LG=1,NPG
  1331.           PDR=PGSQ(LG)*DEUPI*RPG(LG)
  1332.           C1=C1CNG*MELVA1.VELCHE(LG,NK1)*PDR
  1333.          DO     N=1,IDIM
  1334.           I1=I1+1
  1335.           C3=MELVA3.VELCHE((N-1)*NPG+LG,NK3)*C1
  1336.          DO     I=1,MP
  1337.           WSC(I1,I)=WS(I,LG)*C3
  1338.          ENDDO
  1339.          ENDDO
  1340.         ENDDO
  1341.  
  1342.         DO 540 I=1,MP
  1343.         DO 541 J=1,NP
  1344.           U3=0.D0
  1345.         DO 542 I1=1,NIPG
  1346.           U3=U3+WSC(I1,I) * HRD(I1,J)
  1347.  542      CONTINUE
  1348.  
  1349.         DO 546 N=1,NC
  1350.           SM1(I,N)   = SM1(I,N) - U3*TN1(J,N)
  1351.  546    CONTINUE
  1352.  
  1353.  541    CONTINUE
  1354.  540    CONTINUE
  1355.         ENDIF
  1356. C..........Fin Cas CNG uniquement (complément)
  1357.        ENDIF
  1358. C..........Fin Complément pour IKOMP=2 ......................
  1359.  
  1360. C=======================================================================
  1361. c Cas 2D axi pour la QDM
  1362.          IF(IHV.EQ.1.AND.IAXI.EQ.2)THEN
  1363.         DO     LG=1,NPG
  1364.           PDR=PGSQ(LG)*DEUPI/RPG(LG)
  1365.           C2=C2LAPN*MELVA2.VELCHE(LG,NK2)*PDR
  1366. * pv semble inutile et faux         C3=MELVA3.VELCHE((N-1)*NPG+LG,NK3)
  1367.          DO     I=1,MP
  1368.           WTC(LG,I)=WT(I,LG)*C2
  1369.          ENDDO
  1370.         ENDDO
  1371.  
  1372.         DO 310 I=1,MP
  1373.         DO 311 J=1,NP
  1374.          UR=0.D0
  1375.         DO 312 LG=1,NPG
  1376.           UR=UR+WTC(LG,I)*FN(J,LG)
  1377.  312      CONTINUE
  1378.           SM1(I,1) = SM1(I,1) + AI1*UR*TN1(J,1)
  1379.           RF1(J,I,2) = RF1(J,I,1)
  1380.           RF1(J,I,1)=RF1(J,I,1)+(AIEX*UR)
  1381.  311    CONTINUE
  1382.  310    CONTINUE
  1383.          ENDIF
  1384. c Cas 2D axi pour la QDM  Fin
  1385. C=======================================================================
  1386.  
  1387.         ENDIF
  1388.  
  1389. C:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
  1390.         IF(TSOUR) THEN
  1391. C...... Source
  1392.         DO 610 I=1,MP
  1393.          DO 617 N=1,NC
  1394.           U4=0.D0
  1395.           DO 615 LG=1,NPG
  1396.           C4=MELVA4.VELCHE((N-1)*NPG+LG,NK4)
  1397.           U4=U4+WT(I,LG)*PGSQ(LG)*C4*DEUPI*RPG(LG)
  1398.  615      CONTINUE
  1399.           SM1(I,N)=SM1(I,N)+ U4
  1400.  617     CONTINUE
  1401.  610    CONTINUE
  1402. C...... Source Fin
  1403.         ENDIF
  1404. C=======================================================================
  1405.  
  1406. C>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
  1407.       IF(TMATRI.AND.KKMACO.NE.2)THEN
  1408. C ...... Chargement Rigidite ou Matrik
  1409. C Cas RIGIDITE
  1410.           IF(XRIG)THEN
  1411. *           XMATRI=IMATTT(KE)
  1412.            DO     I=1,MP
  1413.            DO     J=1,NP
  1414.            RE(I,J,ke)=RF1(J,I,1)
  1415.            ENDDO
  1416.            ENDDO
  1417. *           SEGDES XMATRI
  1418.           ELSE
  1419. C Cas MATRIK
  1420.            NBMM=1
  1421.            IF(ICAL)NBMM=NBME
  1422.            DO     N=1,NBMM
  1423.            JMATR1=JRIGEL(4,NMATR0+N)
  1424.            IPM4=JMATR1.LIZAFM(L,1)
  1425.            DO     I=1,NP
  1426.            DO     J=1,NP
  1427.            IPM4.AM(KE,J,I)=RF1(J,I,N)
  1428.            ENDDO
  1429.            ENDDO
  1430.            ENDDO
  1431.           ENDIF
  1432.       ENDIF
  1433.  
  1434. C ...... Chargement Second membre
  1435.         DO     I=1,NP
  1436.          I1=LECT(IPT1.NUM(I,KE))
  1437.         DO     N=1,NC
  1438.          MPOVA1.VPOCHA(I1,N)=MPOVA1.VPOCHA(I1,N)+SM1(I,N)
  1439.         ENDDO
  1440.         ENDDO
  1441. C<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
  1442.  
  1443.  108    CONTINUE
  1444.  
  1445.          SEGSUP MCHAM1,MELVA1
  1446.          SEGSUP MCHAM2,MELVA2
  1447.          SEGSUP MCHAM3,MELVA3
  1448.  
  1449.       IF(TSOUR)THEN
  1450.          SEGSUP MCHAM4,MELVA4
  1451.       ENDIF
  1452.  
  1453.       IF(IKOMP.EQ.2)THEN
  1454.          SEGSUP MCHAM5,MELVA5
  1455.       ENDIF
  1456.  
  1457.  
  1458.  
  1459.       SEGDES IPT1,IPT2
  1460.  
  1461.       IF(TMATRI)THEN
  1462. C Cas RIGIDITE
  1463.        IF(XRIG)THEN
  1464.         SEGDES xMATRI
  1465.        ELSE
  1466. C Cas MATRIK
  1467.         SEGDES IPM1
  1468.         IF(NBME.GE.2)SEGDES IPM2
  1469.         IF(NBME.GE.3)SEGDES IPM3
  1470.        ENDIF
  1471.       ENDIF
  1472.  
  1473.       SEGSUP IZFFM,IZHR,IZF1,IZH2
  1474.       SEGSUP SAJT
  1475.  
  1476.  101  CONTINUE
  1477.  
  1478.       IF(TMATRI)THEN
  1479.        IF(.NOT.XRIG)THEN
  1480.         NBMM=JRIGEL(/2)
  1481.         IF(NBMM.NE.0)THEN
  1482.         DO 141 M=1,NBMM
  1483.         JMATRI=JRIGEL(4,M)
  1484.         SEGDES JMATRI
  1485.  141    CONTINUE
  1486.         ENDIF
  1487.        ENDIF
  1488.       ENDIF
  1489.  
  1490. C"""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""
  1491. C"""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""
  1492. C"""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""
  1493.  
  1494.  
  1495.          SEGSUP MCHEL1
  1496.          SEGSUP MCHEL2
  1497.          SEGSUP MCHEL3
  1498.  
  1499.       IF(TSOUR)THEN
  1500.          SEGSUP MCHEL4
  1501.       ENDIF
  1502.  
  1503.       IF(IKOMP.EQ.2)THEN
  1504.          SEGSUP MCHEL5
  1505.       ENDIF
  1506.  
  1507.       SEGDES MCHPO1,MPOVA1
  1508.       SEGDES MELEME,MELEM2
  1509.  
  1510.       SEGSUP MLENTI
  1511. c     IF(TLINCO)THEN
  1512.       SEGSUP MLENT1
  1513.        SEGDES MTETA1
  1514. c     ENDIF
  1515.  
  1516.       IF(TMATRI)THEN
  1517. C Cas RIGIDITE
  1518.        IF(XRIG)THEN
  1519.         SEGDES MRIGID
  1520.        ELSE
  1521. C Cas MATRIK
  1522.         SEGDES MATRIK
  1523.        ENDIF
  1524.       ENDIF
  1525.  
  1526.  
  1527. c     segact MCHPO1
  1528. c     MSOUP1=MCHPO1.IPCHP(1)
  1529. c     segact MSOUP1
  1530. c     write(6,*)MCHPO1.IPCHP(/1),MSOUP1.NOCOMP(/2)
  1531. c     write(6,*)(MSOUP1.NOCOMP(i),i=1,MSOUP1.NOCOMP(/2))
  1532. c     if(ihv.eq.1)call ecmo(mtabz,'TOTO','CHPOINT',MCHPO1)
  1533.  
  1534. C*************************************************************************
  1535.  
  1536.       IF(KKMACO.EQ.1)THEN
  1537.         TYPE='MATRIK'
  1538.         write(6,*)'YCLS ON ECRIT DANS UNE TABLE'
  1539.         CALL ECMO(MTAB1,NMACO,TYPE,MATRIK)
  1540.       ENDIF
  1541.  
  1542.       IF(AIMPL.EQ.0.D0)THEN
  1543.         NRIGE=7
  1544.         NKID =9
  1545.         NKMT =7
  1546.         NMATRI=0
  1547.         SEGINI MATRIK
  1548.       ENDIF
  1549.  
  1550. c      write(6,*)' FIN YCLS'
  1551.       RETURN
  1552.  1001 FORMAT(20(1X,I5))
  1553.  1002 FORMAT(10(1X,1PE11.4))
  1554.       END
  1555.  
  1556.  
  1557.  
  1558.  

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