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    PV090527  26/04/30    21:16:49     12529          
  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.         rigrel=0
  780.         SEGINI xMATRI
  781.         IRIGEL(4,NRI+L)=xMATRI
  782.         xmatri.symre=irigel(7,nri+l)
  783. c       write(6,*)'NELRIG,IMATRI=',NELRIG,IMATRI
  784.  
  785. *        DO 104 K=1,NELRIG
  786. *        SEGINI XMATRI
  787. c       write(6,*)'NLIGRD,NLIGRP,XMATRI=',NLIGRD,NLIGRP,XMATRI
  788. *        IMATTT(K)=XMATRI
  789. * 104    CONTINUE
  790.  
  791.         ELSE
  792. C Cas MATRIK
  793.         SEGINI IPM1
  794.         mtail=(IPM1.AM(/1))*(ipm1.am(/2))*(ipm1.am(/3))
  795.         JMATR1=JRIGEL(4,NMATR0+1)
  796.         JMATR1.LIZAFM(L,1)=IPM1
  797.         IPM2=IPM1
  798.         IPM3=IPM1
  799.          IF(NBME.GE.2)THEN
  800.           JMATR2=JRIGEL(4,NMATR0+2)
  801.           IF(IAXI.NE.0)THEN
  802.            SEGINI IPM2
  803.         mtail=(IPM2.AM(/1))*(ipm2.am(/2))*(ipm2.am(/3))
  804.            JMATR2.LIZAFM(L,1)=IPM2
  805.            ICAL=.TRUE.
  806.           ELSE
  807.            IPM2=IPM1
  808.            JMATR2.LIZAFM(L,1)=IPM2
  809.            ICAL=.FALSE.
  810.           ENDIF
  811.          ENDIF
  812.          IF(NBME.GE.3)THEN
  813.           JMATR3=JRIGEL(4,NMATR0+3)
  814.            IPM3=IPM1
  815.            JMATR3.LIZAFM(L,1)=IPM3
  816.            ICAL=.FALSE.
  817.          ENDIF
  818.        ENDIF
  819.       ENDIF
  820. c       write(6,*)' YCLS appelé par ',NOMPER,' Taille matrice=',mtail
  821.  
  822.       ELSE
  823.  
  824. C Cas RIGIDITE
  825.        IF(XRIG)THEN
  826.  
  827.        ELSE
  828. C Cas MATRIK
  829.         JMATR1=JRIGEL(4,1)
  830.         SEGACT JMATR1
  831.         IPM1=JMATR1.LIZAFM(L,1)
  832.         IPM2=IPM1
  833.         IPM3=IPM1
  834.         SEGACT IPM1
  835.          IF(NBME.GE.2)THEN
  836.           JMATR2=JRIGEL(4,2)
  837.           SEGACT JMATR2
  838.           IPM2=JMATR2.LIZAFM(L,1)
  839.           SEGACT IPM2
  840.          ENDIF
  841.          IF(NBME.GE.3)THEN
  842.           JMATR3=JRIGEL(4,3)
  843.           IPM3=JMATR3.LIZAFM(L,1)
  844.           SEGACT IPM3
  845.          ENDIF
  846.        ENDIF
  847.       ENDIF
  848. C FIN MATRIK . FIN MATRIK . FIN MATRIK . FIN MATRIK . FIN MATRIK . FIN M
  849. C ......................................................................
  850.  
  851.  
  852. C----Ro = 1.
  853.       IK1=1
  854.       MCHAM1=MCHEL1.ICHAML(L)
  855.       SEGACT MCHAM1
  856.  
  857.       MELVA1=MCHAM1.IELVAL(1)
  858.       SEGACT MELVA1
  859.       N1PTEL=MELVA1.VELCHE(/1)
  860.       N1EL=MELVA1.VELCHE(/2)
  861.         IF(N1EL.EQ.1)THEN
  862.         IK1=1
  863.         ELSEIF(N1EL.EQ.NBEL)THEN
  864.         IK1=0
  865.         ENDIF
  866.  
  867. C----Nu ou Lambda
  868.       IK2=1
  869.       MCHAM2=MCHEL2.ICHAML(L)
  870.       SEGACT MCHAM2
  871.  
  872.       MELVA2=MCHAM2.IELVAL(1)
  873.       SEGACT MELVA2
  874.       N1PTEL=MELVA2.VELCHE(/1)
  875.       N1EL=MELVA2.VELCHE(/2)
  876.         IF(N1EL.EQ.1)THEN
  877.         IK2=1
  878.         ELSEIF(N1EL.EQ.NBEL)THEN
  879.         IK2=0
  880.         ENDIF
  881.  
  882.       IF(MCHELG.NE.0)THEN
  883.       MCHAMG=MCHELG.ICHAML(L)
  884.       SEGACT MCHAMG
  885.  
  886.       MELVAG=MCHAMG.IELVAL(1)
  887.       SEGACT MELVAG
  888.       ENDIF
  889.  
  890. C----U
  891.       IK3=1
  892.       MCHAM3=MCHEL3.ICHAML(L)
  893.       SEGACT MCHAM3
  894.  
  895.       MELVA3=MCHAM3.IELVAL(1)
  896.       SEGACT MELVA3
  897.       N1PTEL=MELVA3.VELCHE(/1)
  898.       N1EL=MELVA3.VELCHE(/2)
  899.         IF(N1EL.EQ.1)THEN
  900.         IK3=1
  901.         ELSEIF(N1EL.EQ.NBEL)THEN
  902.         IK3=0
  903.         ENDIF
  904.  
  905. C----sources
  906.       IK4=1
  907.        IF(TSOUR)THEN
  908.       MCHAM4=MCHEL4.ICHAML(L)
  909.       SEGACT MCHAM4
  910.  
  911.       MELVA4=MCHAM4.IELVAL(1)
  912.       SEGACT MELVA4
  913.       N1PTEL=MELVA4.VELCHE(/1)
  914.       N1EL=MELVA4.VELCHE(/2)
  915.         IF(N1EL.EQ.1)THEN
  916.         IK4=1
  917.         ELSEIF(N1EL.EQ.NBEL)THEN
  918.         IK4=0
  919.         ENDIF
  920.        ENDIF
  921.  
  922. C----Div (rocp u )
  923.       IK5=1
  924.        IF(IKOMP.EQ.2)THEN
  925.       MCHAM5=MCHEL5.ICHAML(L)
  926.       SEGACT MCHAM5
  927.  
  928.       MELVA5=MCHAM5.IELVAL(1)
  929.       SEGACT MELVA5
  930.       N1PTEL=MELVA5.VELCHE(/1)
  931.       N1EL=MELVA5.VELCHE(/2)
  932.         IF(N1EL.EQ.1)THEN
  933.         IK5=1
  934.         ELSEIF(N1EL.EQ.NBEL)THEN
  935.         IK5=0
  936.         ENDIF
  937.        ENDIF
  938.  
  939. C      write(6,*)' AVANT 108 NC=',NC,' NBEL=',NBEL,MP,NP,NC
  940. C===============================================
  941.         AI1=AIMPL-1.D0
  942.         IF(AIMPL.EQ.0.D0)THEN
  943.           AI2=-1.D0
  944.         ELSE
  945.           AI2=AI1/AIMPL
  946.         ENDIF
  947.  
  948.         DO 108 KE=1,NBEL
  949.  
  950.         NK1=KE + IK1*(1 - KE)
  951.         NK2=KE + IK2*(1 - KE)
  952.         NK3=KE + IK3*(1 - KE)
  953.         NK4=KE + IK4*(1 - KE)
  954.         NK5=KE + IK5*(1 - KE)
  955.  
  956.         DO     I=1,NP
  957.         J=IPT1.NUM(I,KE)
  958.         DO     N=1,IDIM
  959.         XYZ(N,I)=XCOOR((J-1)*(IDIM+1)+N)
  960.         ENDDO
  961.         ENDDO
  962.  
  963.         DO     I=1,NP
  964.          I1=MLENT1.LECT(IPT1.NUM(I,KE))
  965.         DO     N=1,NC
  966.          TN1(I,N)=MTETA1.VPOCHA(I1,N)
  967.         ENDDO
  968.         ENDDO
  969.  
  970.         CALL CALJBR(FN,GR,PG,XYZ,HR,PGSQ,RPG,NES,
  971.      *  IDIM,NP,NPG,IAXI,AIRE,AJ,SGN)
  972.  
  973.         IF (IDCEN.EQ.0.OR.IDCEN.EQ.1) THEN
  974.         CALL RSETD(WT,FN,(NP*NPG))
  975.         CALL RSETD(WS,FN,(NP*NPG))
  976.         ELSE
  977.         CALL CALDEC(WT,WS,XYZ,GR,HR,FN,NES,IDIM,NP,NPG,AJT,
  978.      &  IDCEN,CMD,MELVA1.VELCHE(1,NK1),MELVA2.VELCHE(1,NK2),
  979.      &  MELVA3.VELCHE(1,NK3),TN1,NC,IKOMP,XREF,AIRE,KE)
  980.         ENDIF
  981.  
  982.         CALL INITD(SM1,(MP*IDIM),0.D0)
  983.  
  984. C:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
  985.         IF(KKMACO.EQ.2)THEN
  986.  
  987. C ...... Chargement Rigidite ou Matrik
  988. C Cas RIGIDITE
  989.           IF(XRIG)THEN
  990. *           XMATRI=IMATTT(KE)
  991.            DO  I=1,MP
  992.            DO  J=1,NP
  993.            RE(I,J,ke)=RF1(J,I,1)
  994.            ENDDO
  995.            ENDDO
  996. *           SEGDES XMATRI
  997.           ELSE
  998. C Cas MATRIK
  999.            DO  N=1,NBME
  1000.            JMATR1=JRIGEL(4,N)
  1001.            IPM4=JMATR1.LIZAFM(L,1)
  1002.            DO  I=1,MP
  1003.            DO  J=1,NP
  1004.            RF1(J,I,N)=IPM4.AM(KE,J,I)
  1005.            ENDDO
  1006.            ENDDO
  1007.            ENDDO
  1008.           ENDIF
  1009. C...... Multiplication Matrice Vecteur
  1010.  
  1011.         DO 710 I=1,MP
  1012.         DO 711 J=1,NP
  1013.  
  1014.         DO 716 N=1,NC
  1015.           SM1(I,N)   = SM1(I,N) + AI2*RF1(J,I,N)*TN1(J,N)
  1016.  716    CONTINUE
  1017.  
  1018.  711    CONTINUE
  1019.  710    CONTINUE
  1020.  
  1021.         ELSEIF(KKMACO.NE.2)THEN
  1022. C.................................................................
  1023. C...... Convection  Laplacien ....  préparation forme conservative
  1024. C...... Convection  Laplacien ....             et non conservative
  1025. C préparation
  1026.           AKOMP=1.D0
  1027.           IF(IKOMP.EQ.1)AKOMP=-1.D0
  1028.           NIPG=NPG*IDIM
  1029.           I1=0
  1030.         DO  LG=1,NPG
  1031.           PDR=PGSQ(LG)*DEUPI*RPG(LG)
  1032.           C1=C1CNG*MELVA1.VELCHE(LG,NK1)*PDR*AKOMP
  1033.           C2=C2LAPN*MELVA2.VELCHE(LG,NK2)*PDR
  1034.          DO  N=1,IDIM
  1035.           I1=I1+1
  1036.           C3=MELVA3.VELCHE((N-1)*NPG+LG,NK3)*C1
  1037.          DO  I=1,MP
  1038.           WTC(I1,I)=WT(I,LG)*C3 + IZWH.HR(N,I,LG)*C2
  1039.          ENDDO
  1040.          ENDDO
  1041.         ENDDO
  1042.  
  1043.           I1=0
  1044.         DO     LG=1,NPG
  1045.          DO     N=1,IDIM
  1046.           I1=I1+1
  1047.          DO     J=1,NP
  1048.           HRD(I1,J)=HR(N,J,LG)
  1049.          ENDDO
  1050.          ENDDO
  1051.         ENDDO
  1052.  
  1053.        IF(IKOMP.EQ.0.OR.IKOMP.EQ.2)THEN
  1054. C...... Convection .... forme non conservative
  1055.         DO     I=1,MP
  1056.         DO     J=1,NP
  1057.           U3=0.D0
  1058.         DO 402 I1=1,NIPG
  1059.           U3=U3+(WTC(I1,I) * HRD(I1,J))
  1060.  402      CONTINUE
  1061.  
  1062.         DO 403 N=1,NC
  1063.           SM1(I,N)   = SM1(I,N) + AI1*U3*TN1(J,N)
  1064.  403    CONTINUE
  1065.           RF1(J,I,1) = AIEX*U3
  1066.  
  1067.         ENDDO
  1068.         ENDDO
  1069.        ENDIF
  1070.  
  1071.        IF(IKOMP.EQ.1)THEN
  1072. C...... Convection .... forme conservative
  1073.         DO     I=1,MP
  1074.         DO     J=1,NP
  1075.           U3=0.D0
  1076.         DO 412 I1=1,NIPG
  1077.           U3=U3+(WTC(I1,I) * HRD(I1,J))
  1078.  412      CONTINUE
  1079.  
  1080.         DO 413 N=1,NC
  1081.           SM1(J,N)   = SM1(J,N) + AI1*U3*TN1(I,N)
  1082.  413    CONTINUE
  1083.           RF1(I,J,1) = AIEX*U3
  1084.  
  1085.         ENDDO
  1086.         ENDDO
  1087.        ENDIF
  1088.  
  1089. C...... Fin Convection ....
  1090.  
  1091. C.......... Cas KMU=1 et FTAU pour QDM (Formulation en To)
  1092.       IF(KMU.EQ.1)THEN
  1093. c       write(6,*)' Cas FTAU'
  1094.  
  1095.         DO      LG=1,NPG
  1096.           PDR=PGSQ(LG)*DEUPI*RPG(LG)
  1097.           C2=C2LAPN*MELVA2.VELCHE(LG,NK2)*PDR
  1098.          DO      N=1,IDIM
  1099.          DO      J=1,NP
  1100.           FNE(LG,J,N)=HR(N,J,LG)*C2
  1101.           HRE(LG,J,N)=HR(N,J,LG)
  1102.          ENDDO
  1103.          ENDDO
  1104.         ENDDO
  1105.  
  1106.       IF(IDIM.EQ.2)THEN
  1107.         DO 4192 I=1,NP
  1108.           US1=0.D0
  1109.           US2=0.D0
  1110.         DO 4182 J=1,NP
  1111.           U11=0.D0
  1112.           U12=0.D0
  1113.           U21=0.D0
  1114.           U22=0.D0
  1115.         DO 4172 LG=1,NPG
  1116.           U11=U11+HRE(LG,J,1)*FNE(LG,I,1)
  1117.           U12=U12+HRE(LG,J,1)*FNE(LG,I,2)
  1118.           U21=U21+HRE(LG,J,2)*FNE(LG,I,1)
  1119.           U22=U22+HRE(LG,J,2)*FNE(LG,I,2)
  1120.  4172   CONTINUE
  1121.           US1 = US1 - U11*TN1(J,1) - U12*TN1(J,2)
  1122.           US2 = US2 - U21*TN1(J,1) - U22*TN1(J,2)
  1123.  4182   CONTINUE
  1124.           SM1(I,1) = SM1(I,1) + US1
  1125.           SM1(I,2) = SM1(I,2) + US2
  1126.  4192   CONTINUE
  1127.       ELSE
  1128.         DO 5193 I=1,NP
  1129.           US1=0.D0
  1130.           US2=0.D0
  1131.           US3=0.D0
  1132.         DO 5183 J=1,NP
  1133.           U11=0.D0
  1134.           U12=0.D0
  1135.           U13=0.D0
  1136.           U21=0.D0
  1137.           U22=0.D0
  1138.           U23=0.D0
  1139.           U31=0.D0
  1140.           U32=0.D0
  1141.           U33=0.D0
  1142.         DO 5173 LG=1,NPG
  1143.           U11=U11+HRE(LG,J,1)*FNE(LG,I,1)
  1144.           U12=U12+HRE(LG,J,1)*FNE(LG,I,2)
  1145.           U13=U13+HRE(LG,J,1)*FNE(LG,I,3)
  1146.  
  1147.           U21=U21+HRE(LG,J,2)*FNE(LG,I,1)
  1148.           U22=U22+HRE(LG,J,2)*FNE(LG,I,2)
  1149.           U23=U23+HRE(LG,J,2)*FNE(LG,I,3)
  1150.  
  1151.           U31=U31+HRE(LG,J,3)*FNE(LG,I,1)
  1152.           U32=U32+HRE(LG,J,3)*FNE(LG,I,2)
  1153.           U33=U33+HRE(LG,J,3)*FNE(LG,I,3)
  1154.  5173   CONTINUE
  1155.           US1 = US1 - U11*TN1(J,1) - U12*TN1(J,2) - U13*TN1(J,3)
  1156.           US2 = US2 - U21*TN1(J,1) - U22*TN1(J,2) - U23*TN1(J,3)
  1157.           US3 = US3 - U31*TN1(J,1) - U32*TN1(J,2) - U33*TN1(J,3)
  1158.  5183   CONTINUE
  1159.           SM1(I,1) = SM1(I,1) + US1
  1160.           SM1(I,2) = SM1(I,2) + US2
  1161.           SM1(I,3) = SM1(I,3) + US3
  1162.  5193   CONTINUE
  1163.       ENDIF
  1164.  
  1165.  
  1166.       ENDIF
  1167. C...... FIN Cas KMU=1 et FTAU pour QDM (Formulation en To)
  1168.  
  1169. C.......... Cas KMU=2 et MUVARI pour QDM (Formulation en Grad de Mu)
  1170.       IF(KMU.EQ.2.AND.MCHELG.NE.0)THEN
  1171. c       write(6,*)' Cas MUVARI'
  1172.  
  1173.           CALL INITD(GMU,(NPG*IDIM),0.D0)
  1174.         DO     LG=1,NPG
  1175.         DO     I=1,NP
  1176.         DO     N=1,NC
  1177.           GMU(LG,N)=GMU(LG,N)+HR(N,I,LG)*MELVAG.VELCHE(I,NK2)
  1178.         ENDDO
  1179.         ENDDO
  1180.         ENDDO
  1181.  
  1182.  
  1183.         DO     LG=1,NPG
  1184.           PDR=PGSQ(LG)*DEUPI*RPG(LG)
  1185.          DO     N=1,IDIM
  1186.          DO     J=1,NP
  1187.           FNE(LG,J,N)=FN(J,LG)*GMU(LG,N)
  1188.           HRE(LG,J,N)=HR(N,J,LG)*PDR
  1189.          ENDDO
  1190.          ENDDO
  1191.          ENDDO
  1192.  
  1193.       IF(IDIM.EQ.2)THEN
  1194.         DO 419 I=1,NP
  1195.           US1=0.D0
  1196.           US2=0.D0
  1197.         DO 418 J=1,NP
  1198.           U11=0.D0
  1199.           U12=0.D0
  1200.           U21=0.D0
  1201.           U22=0.D0
  1202.         DO 417 LG=1,NPG
  1203.           U11=U11+HRE(LG,J,1)*FNE(LG,I,1)
  1204.           U12=U12+HRE(LG,J,1)*FNE(LG,I,2)
  1205.           U21=U21+HRE(LG,J,2)*FNE(LG,I,1)
  1206.           U22=U22+HRE(LG,J,2)*FNE(LG,I,2)
  1207.  417    CONTINUE
  1208.           US1 = US1 - U11*TN1(J,1) - U12*TN1(J,2)
  1209.           US2 = US2 - U21*TN1(J,1) - U22*TN1(J,2)
  1210.  418    CONTINUE
  1211.           SM1(I,1) = SM1(I,1) + US1
  1212.           SM1(I,2) = SM1(I,2) + US2
  1213.  419    CONTINUE
  1214.       ELSE
  1215.         DO 4193 I=1,NP
  1216.           US1=0.D0
  1217.           US2=0.D0
  1218.           US3=0.D0
  1219.         DO 4183 J=1,NP
  1220.           U11=0.D0
  1221.           U12=0.D0
  1222.           U13=0.D0
  1223.           U21=0.D0
  1224.           U22=0.D0
  1225.           U23=0.D0
  1226.           U31=0.D0
  1227.           U32=0.D0
  1228.           U33=0.D0
  1229.         DO 4173 LG=1,NPG
  1230. c         U11=U11+FN(I,LG)*HRE(LG,J,1)*GMU(LG,1)
  1231. c         U12=U12+FN(I,LG)*HRE(LG,J,1)*GMU(LG,2)
  1232. c         U21=U21+FN(I,LG)*HRE(LG,J,2)*GMU(LG,1)
  1233. c         U22=U22+FN(I,LG)*HRE(LG,J,2)*GMU(LG,2)
  1234.  
  1235.           U11=U11+HRE(LG,J,1)*FNE(LG,I,1)
  1236.           U12=U12+HRE(LG,J,1)*FNE(LG,I,2)
  1237.           U13=U13+HRE(LG,J,1)*FNE(LG,I,3)
  1238.  
  1239.           U21=U21+HRE(LG,J,2)*FNE(LG,I,1)
  1240.           U22=U22+HRE(LG,J,2)*FNE(LG,I,2)
  1241.           U23=U23+HRE(LG,J,2)*FNE(LG,I,3)
  1242.  
  1243.           U31=U31+HRE(LG,J,3)*FNE(LG,I,1)
  1244.           U32=U32+HRE(LG,J,3)*FNE(LG,I,2)
  1245.           U33=U33+HRE(LG,J,3)*FNE(LG,I,3)
  1246.  4173   CONTINUE
  1247.           US1 = US1 - U11*TN1(J,1) - U12*TN1(J,2) - U13*TN1(J,3)
  1248.           US2 = US2 - U21*TN1(J,1) - U22*TN1(J,2) - U23*TN1(J,3)
  1249.           US3 = US3 - U31*TN1(J,1) - U32*TN1(J,2) - U33*TN1(J,3)
  1250.  4183   CONTINUE
  1251.           SM1(I,1) = SM1(I,1) + US1
  1252.           SM1(I,2) = SM1(I,2) + US2
  1253.           SM1(I,3) = SM1(I,3) + US3
  1254.  4193   CONTINUE
  1255.       ENDIF
  1256.  
  1257.  
  1258.       ENDIF
  1259. C...... FIN Cas KMU=2 et MUVARI pour QDM (Formulation en Grad de Mu)
  1260.  
  1261. C.......... Cas CNG uniquement (complément)
  1262.         IF(IDCEN.EQ.5)THEN
  1263. C préparation
  1264.           I1=0
  1265.         DO     LG=1,NPG
  1266.           PDR=PGSQ(LG)*DEUPI*RPG(LG)
  1267.           C1=C1CNG*MELVA1.VELCHE(LG,NK1)*PDR
  1268.          DO     N=1,IDIM
  1269.           I1=I1+1
  1270.           C3=MELVA3.VELCHE((N-1)*NPG+LG,NK3)*C1
  1271.          DO     I=1,MP
  1272.           WSC(I1,I)=WS(I,LG)*C3
  1273.          ENDDO
  1274.          ENDDO
  1275.         ENDDO
  1276.  
  1277.         DO 510 I=1,MP
  1278.         DO 511 J=1,NP
  1279.           U3=0.D0
  1280.         DO 512 I1=1,NIPG
  1281.           U3=U3+WSC(I1,I) * HRD(I1,J)
  1282.  512      CONTINUE
  1283.  
  1284.         DO 516 N=1,NC
  1285.           SM1(I,N)   = SM1(I,N) - U3*TN1(J,N)
  1286.  516    CONTINUE
  1287.  
  1288.  511    CONTINUE
  1289.  510    CONTINUE
  1290.         ENDIF
  1291. C..........Fin Cas CNG uniquement (complément)
  1292.  
  1293. C..............................................................
  1294. C..............................................................
  1295.         IF(IKOMP.EQ.2)THEN
  1296. C...... Convection  Laplacien .... Formulation conservative
  1297. C préparation
  1298.         DO     LG=1,NPG
  1299.           PDR=PGSQ(LG)*DEUPI*RPG(LG)
  1300.           C5=C1CNG*MELVA5.VELCHE(LG,NK5)*PDR
  1301.          DO     I=1,MP
  1302.           WTC(LG,I)=WT(I,LG)*C5
  1303.          ENDDO
  1304.         ENDDO
  1305.  
  1306.         DO     LG=1,NPG
  1307.          DO     J=1,NP
  1308.           FND(LG,J)=FN(J,LG)
  1309.          ENDDO
  1310.         ENDDO
  1311.  
  1312.         DO 440 I=1,MP
  1313.         DO 441 J=1,NP
  1314.           U3=0.D0
  1315.         DO 442 LG=1,NPG
  1316.           U3=U3+WTC(LG,I) * FND(LG,J)
  1317.  442      CONTINUE
  1318.  
  1319.         DO 446 N=1,NC
  1320.           SM1(I,N)   = SM1(I,N) + AI1*U3*TN1(J,N)
  1321.  446    CONTINUE
  1322.           RF1(J,I,1)=RF1(J,I,1)+(AIEX*U3)
  1323.  
  1324.  441    CONTINUE
  1325.  440    CONTINUE
  1326.  
  1327. C.......... Cas CNG uniquement (complément)
  1328.         IF(IDCEN.EQ.5)THEN
  1329. C préparation
  1330.           I1=0
  1331.         DO     LG=1,NPG
  1332.           PDR=PGSQ(LG)*DEUPI*RPG(LG)
  1333.           C1=C1CNG*MELVA1.VELCHE(LG,NK1)*PDR
  1334.          DO     N=1,IDIM
  1335.           I1=I1+1
  1336.           C3=MELVA3.VELCHE((N-1)*NPG+LG,NK3)*C1
  1337.          DO     I=1,MP
  1338.           WSC(I1,I)=WS(I,LG)*C3
  1339.          ENDDO
  1340.          ENDDO
  1341.         ENDDO
  1342.  
  1343.         DO 540 I=1,MP
  1344.         DO 541 J=1,NP
  1345.           U3=0.D0
  1346.         DO 542 I1=1,NIPG
  1347.           U3=U3+WSC(I1,I) * HRD(I1,J)
  1348.  542      CONTINUE
  1349.  
  1350.         DO 546 N=1,NC
  1351.           SM1(I,N)   = SM1(I,N) - U3*TN1(J,N)
  1352.  546    CONTINUE
  1353.  
  1354.  541    CONTINUE
  1355.  540    CONTINUE
  1356.         ENDIF
  1357. C..........Fin Cas CNG uniquement (complément)
  1358.        ENDIF
  1359. C..........Fin Complément pour IKOMP=2 ......................
  1360.  
  1361. C=======================================================================
  1362. c Cas 2D axi pour la QDM
  1363.          IF(IHV.EQ.1.AND.IAXI.EQ.2)THEN
  1364.         DO     LG=1,NPG
  1365.           PDR=PGSQ(LG)*DEUPI/RPG(LG)
  1366.           C2=C2LAPN*MELVA2.VELCHE(LG,NK2)*PDR
  1367. * pv semble inutile et faux         C3=MELVA3.VELCHE((N-1)*NPG+LG,NK3)
  1368.          DO     I=1,MP
  1369.           WTC(LG,I)=WT(I,LG)*C2
  1370.          ENDDO
  1371.         ENDDO
  1372.  
  1373.         DO 310 I=1,MP
  1374.         DO 311 J=1,NP
  1375.          UR=0.D0
  1376.         DO 312 LG=1,NPG
  1377.           UR=UR+WTC(LG,I)*FN(J,LG)
  1378.  312      CONTINUE
  1379.           SM1(I,1) = SM1(I,1) + AI1*UR*TN1(J,1)
  1380.           RF1(J,I,2) = RF1(J,I,1)
  1381.           RF1(J,I,1)=RF1(J,I,1)+(AIEX*UR)
  1382.  311    CONTINUE
  1383.  310    CONTINUE
  1384.          ENDIF
  1385. c Cas 2D axi pour la QDM  Fin
  1386. C=======================================================================
  1387.  
  1388.         ENDIF
  1389.  
  1390. C:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
  1391.         IF(TSOUR) THEN
  1392. C...... Source
  1393.         DO 610 I=1,MP
  1394.          DO 617 N=1,NC
  1395.           U4=0.D0
  1396.           DO 615 LG=1,NPG
  1397.           C4=MELVA4.VELCHE((N-1)*NPG+LG,NK4)
  1398.           U4=U4+WT(I,LG)*PGSQ(LG)*C4*DEUPI*RPG(LG)
  1399.  615      CONTINUE
  1400.           SM1(I,N)=SM1(I,N)+ U4
  1401.  617     CONTINUE
  1402.  610    CONTINUE
  1403. C...... Source Fin
  1404.         ENDIF
  1405. C=======================================================================
  1406.  
  1407. C>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
  1408.       IF(TMATRI.AND.KKMACO.NE.2)THEN
  1409. C ...... Chargement Rigidite ou Matrik
  1410. C Cas RIGIDITE
  1411.           IF(XRIG)THEN
  1412. *           XMATRI=IMATTT(KE)
  1413.            DO     I=1,MP
  1414.            DO     J=1,NP
  1415.            RE(I,J,ke)=RF1(J,I,1)
  1416.            ENDDO
  1417.            ENDDO
  1418. *           SEGDES XMATRI
  1419.           ELSE
  1420. C Cas MATRIK
  1421.            NBMM=1
  1422.            IF(ICAL)NBMM=NBME
  1423.            DO     N=1,NBMM
  1424.            JMATR1=JRIGEL(4,NMATR0+N)
  1425.            IPM4=JMATR1.LIZAFM(L,1)
  1426.            DO     I=1,NP
  1427.            DO     J=1,NP
  1428.            IPM4.AM(KE,J,I)=RF1(J,I,N)
  1429.            ENDDO
  1430.            ENDDO
  1431.            ENDDO
  1432.           ENDIF
  1433.       ENDIF
  1434.  
  1435. C ...... Chargement Second membre
  1436.         DO     I=1,NP
  1437.          I1=LECT(IPT1.NUM(I,KE))
  1438.         DO     N=1,NC
  1439.          MPOVA1.VPOCHA(I1,N)=MPOVA1.VPOCHA(I1,N)+SM1(I,N)
  1440.         ENDDO
  1441.         ENDDO
  1442. C<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
  1443.  
  1444.  108    CONTINUE
  1445.  
  1446.          SEGSUP MCHAM1,MELVA1
  1447.          SEGSUP MCHAM2,MELVA2
  1448.          SEGSUP MCHAM3,MELVA3
  1449.  
  1450.       IF(TSOUR)THEN
  1451.          SEGSUP MCHAM4,MELVA4
  1452.       ENDIF
  1453.  
  1454.       IF(IKOMP.EQ.2)THEN
  1455.          SEGSUP MCHAM5,MELVA5
  1456.       ENDIF
  1457.  
  1458.  
  1459.  
  1460.       SEGDES IPT1,IPT2
  1461.  
  1462.       IF(TMATRI)THEN
  1463. C Cas RIGIDITE
  1464.        IF(XRIG)THEN
  1465.         SEGDES xMATRI
  1466.        ELSE
  1467. C Cas MATRIK
  1468.         SEGDES IPM1
  1469.         IF(NBME.GE.2)SEGDES IPM2
  1470.         IF(NBME.GE.3)SEGDES IPM3
  1471.        ENDIF
  1472.       ENDIF
  1473.  
  1474.       SEGSUP IZFFM,IZHR,IZF1,IZH2
  1475.       SEGSUP SAJT
  1476.  
  1477.  101  CONTINUE
  1478.  
  1479.       IF(TMATRI)THEN
  1480.        IF(.NOT.XRIG)THEN
  1481.         NBMM=JRIGEL(/2)
  1482.         IF(NBMM.NE.0)THEN
  1483.         DO 141 M=1,NBMM
  1484.         JMATRI=JRIGEL(4,M)
  1485.         SEGDES JMATRI
  1486.  141    CONTINUE
  1487.         ENDIF
  1488.        ENDIF
  1489.       ENDIF
  1490.  
  1491. C"""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""
  1492. C"""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""
  1493. C"""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""
  1494.  
  1495.  
  1496.          SEGSUP MCHEL1
  1497.          SEGSUP MCHEL2
  1498.          SEGSUP MCHEL3
  1499.  
  1500.       IF(TSOUR)THEN
  1501.          SEGSUP MCHEL4
  1502.       ENDIF
  1503.  
  1504.       IF(IKOMP.EQ.2)THEN
  1505.          SEGSUP MCHEL5
  1506.       ENDIF
  1507.  
  1508.       SEGDES MCHPO1,MPOVA1
  1509.       SEGDES MELEME,MELEM2
  1510.  
  1511.       SEGSUP MLENTI
  1512. c     IF(TLINCO)THEN
  1513.       SEGSUP MLENT1
  1514.        SEGDES MTETA1
  1515. c     ENDIF
  1516.  
  1517.       IF(TMATRI)THEN
  1518. C Cas RIGIDITE
  1519.        IF(XRIG)THEN
  1520.         SEGDES MRIGID
  1521.        ELSE
  1522. C Cas MATRIK
  1523.         SEGDES MATRIK
  1524.        ENDIF
  1525.       ENDIF
  1526.  
  1527.  
  1528. c     segact MCHPO1
  1529. c     MSOUP1=MCHPO1.IPCHP(1)
  1530. c     segact MSOUP1
  1531. c     write(6,*)MCHPO1.IPCHP(/1),MSOUP1.NOCOMP(/2)
  1532. c     write(6,*)(MSOUP1.NOCOMP(i),i=1,MSOUP1.NOCOMP(/2))
  1533. c     if(ihv.eq.1)call ecmo(mtabz,'TOTO','CHPOINT',MCHPO1)
  1534.  
  1535. C*************************************************************************
  1536.  
  1537.       IF(KKMACO.EQ.1)THEN
  1538.         TYPE='MATRIK'
  1539.         write(6,*)'YCLS ON ECRIT DANS UNE TABLE'
  1540.         CALL ECMO(MTAB1,NMACO,TYPE,MATRIK)
  1541.       ENDIF
  1542.  
  1543.       IF(AIMPL.EQ.0.D0)THEN
  1544.         NRIGE=7
  1545.         NKID =9
  1546.         NKMT =7
  1547.         NMATRI=0
  1548.         SEGINI MATRIK
  1549.       ENDIF
  1550.  
  1551. c      write(6,*)' FIN YCLS'
  1552.       RETURN
  1553.  1001 FORMAT(20(1X,I5))
  1554.  1002 FORMAT(10(1X,1PE11.4))
  1555.       END
  1556.  
  1557.  
  1558.  
  1559.  
  1560.  
  1561.  

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