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