Télécharger la plaquette Cast3M


fr en

Télécharger chist.eso

Retour à la liste

Numérotation des lignes : 

chist

  1. C CHIST     SOURCE    PV        17/12/08    21:15:42     9660           
  2.       SUBROUTINE CHIST(ipnua1,ipnua2,wrk52,motcas)
  3. *---------------------------------------
  4. * creer un nuage d interpretation du TRC
  5. * pour le modele CEREM (Martinez 1999)
  6. * creer un tableau pour le modele de Leblond (chauffage)
  7. *----------------------------------------
  8.       IMPLICIT INTEGER(I-N)
  9.       IMPLICIT REAL*8 (A-H,O-Z)
  10. -INC PPARAM
  11. -INC CCOPTIO
  12. -INC SMNUAGE
  13. -INC SMLENTI
  14. -INC SMLREEL
  15. -INC DECHE
  16.  
  17.       character*16 motcas
  18.       parameter(nmcer=7,nmtrc=10,nmcha=3)
  19.       character*4 motcer(nmcer),moh,mottrc(nmtrc)
  20.       character*4 motcha(nmcha),motmon(nmcha)
  21.       CHARACTER  TRC1*10
  22.       DIMENSION TE(50),CK(50),CL(50)
  23.       integer indcer(nmcer)
  24.       data motcer/'TP','ZFF','ZFB','TDF','TFF','TDB','TFB'/
  25.       data mottrc/'T   ','TP  ','ZA  ','ZF  ','ZB  ','ZM  ','ZP  ',
  26.      &'ZFP ','ZBP ','MS '/
  27.       data motcha/'TE  ','CK  ','CL  '/
  28.       data motmon/'T   ','CKA ','CLA '/
  29.  
  30.       MNUAG1 = ipnua1
  31.       interr(1) = ipnua1
  32.       if (ipnua1.le.0.and.motcas.eq.'CEREMREFR') then
  33.         moterr(1:8) = 'NHTR'
  34.         call erreur(679)
  35.         return
  36.       endif
  37.       if (ipnua1.le.0.and.motcas.eq.'CEREMCHAU') then
  38.         call erreur(679)
  39.         moterr(1:8) = 'NLEB'
  40.         return
  41.       endif
  42.       segact mnuag1*nomod
  43.       nvar1 = mnuag1.nuanom(/2)
  44.  
  45.       AC1 = valma0(1)
  46.       AR1=valma0(2)
  47.       TMS0=valma0(3)
  48.       BETA=valma0(4)
  49.       AC=valma0(5)
  50.       Aa=valma0(6)
  51.       ZS=valma0(7)
  52.       TPLM=valma0(8)
  53.       C0=valma0(9)
  54.       ALPHA=valma0(10)
  55.       DG0=valma0(11)
  56.       AG=valma0(12)
  57.             ti = valma0(13)
  58.             tf = valma0(14)
  59.             dtemp = valma0(15)
  60.  
  61.       if (motcas.eq.'CEREMREFR') goto 10
  62.       if (motcas.eq.'CEREMCHAU') goto 20
  63. * motcas non prevu
  64.        goto 900
  65.  
  66.  10   continue
  67. * rechercher les composantes, remplir le tableau d adresses
  68.       do im = 1, nmcer
  69.         moh = motcer(im)
  70.         moterr(1:4) = moh
  71.             IMOT=0
  72.             DO 11 IA=1,nvar1
  73.             if (MOH.EQ.mnuag1.nuanom(IA))   then
  74.               imot = ia
  75.               goto 12
  76.             endif
  77.    11       CONTINUE
  78.    12       continue
  79.         if (imot.eq.0) goto 900
  80.         if (( mnuag1.nuatyp(imot).ne.'REAL*8')
  81.      &.AND.( mnuag1.nuatyp(imot).ne.'FLOTTANT')) goto 900
  82.         indcer(im) = mnuag1.nuapoi(imot)
  83.       enddo
  84.       do im =1,nmcer
  85.         nuavf1 = indcer(im)
  86.         segact nuavf1*nomod
  87.       enddo
  88.       segdes mnuag1
  89.       nhist = nuavf1.nuaflo(/1)
  90.  
  91. * simule une table avec des listentiers et des listreels
  92.       jg = nhist
  93.       segini mlenti
  94.       ipnua2 = mlenti
  95.  
  96. C dimensionne les segments
  97.       nbcoup = int ((ti - tf)/dtemp) + 1
  98.  
  99. * traiter chaque champ de donnees
  100.       do ihist = 1,nhist
  101.  
  102. * composantes du nuage initial TP,ZFF,ZFB,TDF,TFF,TDB,TFB
  103.         nuavf1 = indcer(1)
  104.         TP = nuavf1.nuaflo(ihist)
  105.         nuavf1 = indcer(2)
  106.         ZFF =  nuavf1.nuaflo(ihist)
  107.         nuavf1 = indcer(3)
  108.         ZFB =  nuavf1.nuaflo(ihist)
  109.         nuavf1 = indcer(4)
  110.         TDF =  nuavf1.nuaflo(ihist)
  111.         nuavf1 = indcer(5)
  112.         TFF =  nuavf1.nuaflo(ihist)
  113.         nuavf1 = indcer(6)
  114.         TDB =  nuavf1.nuaflo(ihist)
  115.         nuavf1 = indcer(7)
  116.         TFB =  nuavf1.nuaflo(ihist)
  117.  
  118. * on va ranger dans un listenti
  119.         jg = nbcoup + 1
  120.         segini mlent1
  121.       lect(ihist) = mlent1
  122. C
  123. C  INITIALISATION
  124. C
  125.       T = DBLE(TI)
  126.       Z = 1.D0
  127.       ZF = 0.D0
  128.       ZB = 0.D0
  129.       ZM = 0.D0
  130.       ZP = 0.D0
  131.       ZFP = 0.D0
  132.       ZBP = 0.D0
  133. C PARAM POUR LA SAISIE DES COURBES
  134.       DTF = TDF - TFF
  135.       DTB = TDB - TFB
  136.       IUBF1 = IUP (DTB,DTF)
  137.       IUFB1 = IUPS(DTF,DTB)
  138.       DT = DTF/10.D0*IUBF1 + DTB/10.D0*IUFB1
  139.  
  140.  
  141.       icoup = 0
  142. C
  143. c      DO 2101 WHILE (T.GE.DBLE(TF))
  144. C
  145. 2102  CONTINUE
  146.       IF (T.GE.DBLE(TF)) THEN
  147.        icoup = icoup + 1
  148.       IUF1S = IUPS(T, (TFF+DT))
  149.       IUF1  = IUP((TDF-DT), T)
  150.       IUF2S = IUPS(T, (TFF-DT))
  151.       IUF2  = IUP((TFF+DT), T)
  152.       IUF3S = IUPS(T, (TDF-DT))
  153.       IUF3  = IUP((TDF+DT), T)
  154.       IUF4 = IUP((TFF-DT), T)
  155.       ZF = ZFF*(T - TDF)/(TFF -TDF)*IUF1S*IUF1
  156.      & + ZFF*(DT/(TFF-TDF)*(T-(TFF+DT))/(2*DT)+(1+DT/(TFF-TDF)))
  157.      & *IUF2S*IUF2
  158.      & + ZFF*DT/(TFF-TDF)*(T-(TDF+DT))/(2*DT)*IUF3S*IUF3
  159.      & + ZFF*IUF4
  160. C
  161.       IUB1S = IUPS(T, (TFB+DT))
  162.       IUB1  = IUP((TDB-DT), T)
  163.       IUB2S = IUPS(T, (TFB-DT))
  164.       IUB2  = IUP((TFB+DT), T)
  165.       IUB3S = IUPS(T, (TDB-DT))
  166.       IUB3  = IUP((TDB+DT), T)
  167.       IUB4 = IUP((TFB-DT), T)
  168.       ZB = ZFB*(T - TDB)/(TFB -TDB)*IUB1S*IUB1
  169.      & + ZFB*(DT/(TFB-TDB)*(T-(TFB+DT))/(2*DT)+(1+DT/(TFB-TDB)))
  170.      & *IUB2S*IUB2
  171.      & + ZFB*DT/(TFB-TDB)*(T-(TDB+DT))/(2*DT)*IUB3S*IUB3
  172.      & + ZFB*IUB4
  173. C
  174.       ZFP = ZFF/(TFF-TDF)*TP*(IUF1S*IUF1+0.5*(IUF2S*IUF2+IUF3S*IUF3))
  175. C
  176.       ZBP = ZFB/(TFB-TDB)*TP*(IUB1S*IUB1+0.5*(IUB2S*IUB2+IUB3S*IUB3))
  177. C
  178.       Z1 = ZF + ZB
  179.       Z2 = POSITI(Z1, ZS)
  180. * plg
  181.       A = 0
  182. * plg
  183.       TMS = TMS0 + A*Z2
  184. C
  185. C      IF ((TP.LE.TPLM).AND.(T.LE.TFB)) THEN
  186. C      T1 = POSITI(TMS, T)
  187. C      ZM = (1.D0 - ZF - ZB)*(1.D0 - EXP(BETA*T1))
  188. C      ENDIF
  189. C
  190. C      IF ((TP.LE.TPLM).AND.(T.LE.TFB)) THEN
  191. C      IU1 = IUP(TMS, T)
  192. C      ZP = -1.D0*(ZFP + ZBP
  193. C     &    + BETA*(1.D0 - ZF - ZB)*EXP(BETA*T1)*TP*IU1)
  194. C      ELSE
  195.       ZP = -1.D0*(ZFP + ZBP)
  196. C      ENDIF
  197. C
  198.       Z = 1.D0 - ZF - ZB
  199. C
  200. *      WRITE(1,2013) ZFP, ZBP, TMS
  201.          jg = nmtrc
  202.         segini mlreel
  203.        prog(1) = T
  204.        prog(2) = TP
  205.        prog(3) = Z
  206.        prog(4) = ZF
  207.        prog(5) = ZB
  208.        prog(7) = ZP
  209.        prog(8) = ZFP
  210.        prog(9) = ZBP
  211.        prog(10) = TMS
  212.         segdes mlreel
  213.         mlent1.lect(1+icoup) = mlreel
  214.  
  215.  2012 FORMAT(F8.2,1X,F8.2,1X,F9.6,1X,F9.6,1X,F9.6,1X,F9.6)
  216.  2013 FORMAT(F9.6,1X,F9.6,1X,F8.2)
  217. C
  218.       T = T - DBLE(dtemp)
  219. C
  220. c2101  CONTINUE
  221.       GOTO 2102
  222.       END IF
  223. C
  224. c      WRITE(*,*) 'HISTOIRE ENREGISTREE', icoup
  225. C
  226.        segdes mlent1
  227.  3001 CONTINUE
  228.  
  229.       enddo
  230.  
  231.       goto 1000
  232.  
  233. * cas du chauffage
  234.  20   continue
  235. * rechercher les composantes, remplir le tableau d adresses
  236.       do im = 1, nmcha
  237.         moh = motcha(im)
  238.         moterr(1:4) = moh
  239.         call PLACE(mnuag1.nuanom,nvar1,IMOT,MOH)
  240.         if (imot.eq.0) goto 900
  241.         if ( (mnuag1.nuatyp(imot).ne.'REAL*8')
  242.      &.AND.( mnuag1.nuatyp(imot).ne.'FLOTTANT')) goto 900
  243.         indcer(im) = mnuag1.nuapoi(imot)
  244.       enddo
  245.       do im =1,nmcha
  246.         nuavf1 = indcer(im)
  247.         segact nuavf1*nomod
  248.       enddo
  249.       segdes mnuag1
  250.       N = nuavf1.nuaflo(/1)
  251.  
  252. C dimensionne le listenti
  253.       nbcoup = int ((ti - tf)/dtemp) + 1
  254.  
  255.       jg = nbcoup
  256.       segini mlenti
  257.       ipnua2 = mlenti
  258.  
  259. * prepare les tableaux de calcul
  260. C      N=0
  261. C      DO 4001 I=1,50
  262.       DO 4001 I=1,N
  263. c        READ(2,'(f8.2,2(1x,f8.4))',END=4001) TE(I),CK(I),CL(I)
  264. *        READ(2,*,END=4001) TE(I),CK(I),CL(I)
  265.         nuavf1 = indcer(1)
  266.         TE(I) = nuavf1.nuaflo(I)
  267.         nuavf1 = indcer(2)
  268.         CK(I) = nuavf1.nuaflo(I)
  269.         nuavf1 = indcer(3)
  270.         CL(I) = nuavf1.nuaflo(I)
  271. C        N=N+1
  272. 4001  CONTINUE
  273.  
  274. * initialise
  275.       T=DBLE(TI)
  276.  
  277. * traiter chaque champ de donnees
  278.       do icoup = 1,nbcoup
  279.  
  280. c      DO 4011 WHILE (T.GE.DBLE(TF))
  281. 4009  CONTINUE
  282.       IF (T.GE.DBLE(TF)) THEN
  283.         CKA=CK(1)*IUPS(TE(1),T)
  284.         CLA=CL(1)*IUPS(TE(1),T)
  285.         DO 4012 I=1,(N-1)
  286.           CKA=CKA
  287.      &        +((CK(I)*(TE(I+1)-T)+CK(I+1)*(T-TE(I)))/(TE(I+1)-TE(I)))
  288.      &        *IUPS(TE(I+1),T)*IUP(T,TE(I))
  289.           CLA=CLA
  290.      &        +((CL(I)*(TE(I+1)-T)+CL(I+1)*(T-TE(I)))/(TE(I+1)-TE(I)))
  291.      &        *IUPS(TE(I+1),T)*IUP(T,TE(I))
  292. 4012    CONTINUE
  293.         CKA=CKA+CK(N)*IUP(T,TE(N))
  294.         CLA=CLA+CL(N)*IUP(T,TE(N))
  295. *        WRITE(1,'(f8.2,2(1x,f8.4))') T,CKA,CLA
  296.  
  297.        jg =nmcha
  298.        segini mlreel
  299.        lect(icoup) = mlreel
  300.          prog(1) = T
  301.          prog(2) = CKA
  302.          prog(3) = CLA
  303.         segdes mlreel
  304.  
  305.         T=T-DBLE(dtemp)
  306. c4011  CONTINUE
  307. *      GOTO 4009
  308.       END IF
  309.  
  310.       enddo
  311.  
  312.        segdes mlenti
  313.  
  314. c      WRITE(*,*)'CHAUFFAGE ENREGISTRE', t,tf
  315.       goto 1000
  316.  
  317.  
  318.  900  continue
  319.       segdes mnuag1
  320. * la composante %m1:4 du nuage %i1 est absente ou de type incorrect
  321.       call erreur(921)
  322.       RETURN
  323.  
  324.  1000 continue
  325.       do ivar = 1,nvar1
  326.         nuavf1 = indcer(ivar)
  327.         segdes nuavf1
  328.       enddo
  329.       return
  330.  
  331.       END
  332.  
  333.  
  334.  
  335.  
  336.  
  337.  
  338.  
  339.  
  340.  
  341.  
  342.  
  343.  

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