Télécharger chist.eso

Retour à la liste

Numérotation des lignes :

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

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