C CHIST SOURCE PV 17/12/08 21:15:42 9660 SUBROUTINE CHIST(ipnua1,ipnua2,wrk52,motcas) *--------------------------------------- * creer un nuage d interpretation du TRC * pour le modele CEREM (Martinez 1999) * creer un tableau pour le modele de Leblond (chauffage) *---------------------------------------- IMPLICIT INTEGER(I-N) IMPLICIT REAL*8 (A-H,O-Z) -INC PPARAM -INC CCOPTIO -INC SMNUAGE -INC SMLENTI -INC SMLREEL -INC DECHE character*16 motcas parameter(nmcer=7,nmtrc=10,nmcha=3) character*4 motcer(nmcer),moh,mottrc(nmtrc) character*4 motcha(nmcha),motmon(nmcha) CHARACTER TRC1*10 DIMENSION TE(50),CK(50),CL(50) integer indcer(nmcer) data motcer/'TP','ZFF','ZFB','TDF','TFF','TDB','TFB'/ data mottrc/'T ','TP ','ZA ','ZF ','ZB ','ZM ','ZP ', &'ZFP ','ZBP ','MS '/ data motcha/'TE ','CK ','CL '/ data motmon/'T ','CKA ','CLA '/ MNUAG1 = ipnua1 interr(1) = ipnua1 if (ipnua1.le.0.and.motcas.eq.'CEREMREFR') then moterr(1:8) = 'NHTR' call erreur(679) return endif if (ipnua1.le.0.and.motcas.eq.'CEREMCHAU') then call erreur(679) moterr(1:8) = 'NLEB' return endif segact mnuag1*nomod nvar1 = mnuag1.nuanom(/2) AC1 = valma0(1) AR1=valma0(2) TMS0=valma0(3) BETA=valma0(4) AC=valma0(5) Aa=valma0(6) ZS=valma0(7) TPLM=valma0(8) C0=valma0(9) ALPHA=valma0(10) DG0=valma0(11) AG=valma0(12) ti = valma0(13) tf = valma0(14) dtemp = valma0(15) if (motcas.eq.'CEREMREFR') goto 10 if (motcas.eq.'CEREMCHAU') goto 20 * motcas non prevu goto 900 10 continue * rechercher les composantes, remplir le tableau d adresses do im = 1, nmcer moh = motcer(im) moterr(1:4) = moh IMOT=0 DO 11 IA=1,nvar1 if (MOH.EQ.mnuag1.nuanom(IA)) then imot = ia goto 12 endif 11 CONTINUE 12 continue if (imot.eq.0) goto 900 if (( mnuag1.nuatyp(imot).ne.'REAL*8') &.AND.( mnuag1.nuatyp(imot).ne.'FLOTTANT')) goto 900 indcer(im) = mnuag1.nuapoi(imot) enddo do im =1,nmcer nuavf1 = indcer(im) segact nuavf1*nomod enddo segdes mnuag1 nhist = nuavf1.nuaflo(/1) * simule une table avec des listentiers et des listreels jg = nhist segini mlenti ipnua2 = mlenti C dimensionne les segments nbcoup = int ((ti - tf)/dtemp) + 1 * traiter chaque champ de donnees do ihist = 1,nhist * composantes du nuage initial TP,ZFF,ZFB,TDF,TFF,TDB,TFB nuavf1 = indcer(1) TP = nuavf1.nuaflo(ihist) nuavf1 = indcer(2) ZFF = nuavf1.nuaflo(ihist) nuavf1 = indcer(3) ZFB = nuavf1.nuaflo(ihist) nuavf1 = indcer(4) TDF = nuavf1.nuaflo(ihist) nuavf1 = indcer(5) TFF = nuavf1.nuaflo(ihist) nuavf1 = indcer(6) TDB = nuavf1.nuaflo(ihist) nuavf1 = indcer(7) TFB = nuavf1.nuaflo(ihist) * on va ranger dans un listenti jg = nbcoup + 1 segini mlent1 lect(ihist) = mlent1 C C INITIALISATION C T = DBLE(TI) Z = 1.D0 ZF = 0.D0 ZB = 0.D0 ZM = 0.D0 ZP = 0.D0 ZFP = 0.D0 ZBP = 0.D0 C PARAM POUR LA SAISIE DES COURBES DTF = TDF - TFF DTB = TDB - TFB IUBF1 = IUP (DTB,DTF) IUFB1 = IUPS(DTF,DTB) DT = DTF/10.D0*IUBF1 + DTB/10.D0*IUFB1 icoup = 0 C c DO 2101 WHILE (T.GE.DBLE(TF)) C 2102 CONTINUE IF (T.GE.DBLE(TF)) THEN icoup = icoup + 1 IUF1S = IUPS(T, (TFF+DT)) IUF1 = IUP((TDF-DT), T) IUF2S = IUPS(T, (TFF-DT)) IUF2 = IUP((TFF+DT), T) IUF3S = IUPS(T, (TDF-DT)) IUF3 = IUP((TDF+DT), T) IUF4 = IUP((TFF-DT), T) ZF = ZFF*(T - TDF)/(TFF -TDF)*IUF1S*IUF1 & + ZFF*(DT/(TFF-TDF)*(T-(TFF+DT))/(2*DT)+(1+DT/(TFF-TDF))) & *IUF2S*IUF2 & + ZFF*DT/(TFF-TDF)*(T-(TDF+DT))/(2*DT)*IUF3S*IUF3 & + ZFF*IUF4 C IUB1S = IUPS(T, (TFB+DT)) IUB1 = IUP((TDB-DT), T) IUB2S = IUPS(T, (TFB-DT)) IUB2 = IUP((TFB+DT), T) IUB3S = IUPS(T, (TDB-DT)) IUB3 = IUP((TDB+DT), T) IUB4 = IUP((TFB-DT), T) ZB = ZFB*(T - TDB)/(TFB -TDB)*IUB1S*IUB1 & + ZFB*(DT/(TFB-TDB)*(T-(TFB+DT))/(2*DT)+(1+DT/(TFB-TDB))) & *IUB2S*IUB2 & + ZFB*DT/(TFB-TDB)*(T-(TDB+DT))/(2*DT)*IUB3S*IUB3 & + ZFB*IUB4 C ZFP = ZFF/(TFF-TDF)*TP*(IUF1S*IUF1+0.5*(IUF2S*IUF2+IUF3S*IUF3)) C ZBP = ZFB/(TFB-TDB)*TP*(IUB1S*IUB1+0.5*(IUB2S*IUB2+IUB3S*IUB3)) C Z1 = ZF + ZB Z2 = POSITI(Z1, ZS) * plg A = 0 * plg TMS = TMS0 + A*Z2 C C IF ((TP.LE.TPLM).AND.(T.LE.TFB)) THEN C T1 = POSITI(TMS, T) C ZM = (1.D0 - ZF - ZB)*(1.D0 - EXP(BETA*T1)) C ENDIF C C IF ((TP.LE.TPLM).AND.(T.LE.TFB)) THEN C IU1 = IUP(TMS, T) C ZP = -1.D0*(ZFP + ZBP C & + BETA*(1.D0 - ZF - ZB)*EXP(BETA*T1)*TP*IU1) C ELSE ZP = -1.D0*(ZFP + ZBP) C ENDIF C Z = 1.D0 - ZF - ZB C * WRITE(1,2013) ZFP, ZBP, TMS jg = nmtrc segini mlreel prog(1) = T prog(2) = TP prog(3) = Z prog(4) = ZF prog(5) = ZB prog(7) = ZP prog(8) = ZFP prog(9) = ZBP prog(10) = TMS segdes mlreel mlent1.lect(1+icoup) = mlreel 2012 FORMAT(F8.2,1X,F8.2,1X,F9.6,1X,F9.6,1X,F9.6,1X,F9.6) 2013 FORMAT(F9.6,1X,F9.6,1X,F8.2) C T = T - DBLE(dtemp) C c2101 CONTINUE GOTO 2102 END IF C c WRITE(*,*) 'HISTOIRE ENREGISTREE', icoup C segdes mlent1 3001 CONTINUE enddo goto 1000 * cas du chauffage 20 continue * rechercher les composantes, remplir le tableau d adresses do im = 1, nmcha moh = motcha(im) moterr(1:4) = moh call PLACE(mnuag1.nuanom,nvar1,IMOT,MOH) if (imot.eq.0) goto 900 if ( (mnuag1.nuatyp(imot).ne.'REAL*8') &.AND.( mnuag1.nuatyp(imot).ne.'FLOTTANT')) goto 900 indcer(im) = mnuag1.nuapoi(imot) enddo do im =1,nmcha nuavf1 = indcer(im) segact nuavf1*nomod enddo segdes mnuag1 N = nuavf1.nuaflo(/1) C dimensionne le listenti nbcoup = int ((ti - tf)/dtemp) + 1 jg = nbcoup segini mlenti ipnua2 = mlenti * prepare les tableaux de calcul C N=0 C DO 4001 I=1,50 DO 4001 I=1,N c READ(2,'(f8.2,2(1x,f8.4))',END=4001) TE(I),CK(I),CL(I) * READ(2,*,END=4001) TE(I),CK(I),CL(I) nuavf1 = indcer(1) TE(I) = nuavf1.nuaflo(I) nuavf1 = indcer(2) CK(I) = nuavf1.nuaflo(I) nuavf1 = indcer(3) CL(I) = nuavf1.nuaflo(I) C N=N+1 4001 CONTINUE * initialise T=DBLE(TI) * traiter chaque champ de donnees do icoup = 1,nbcoup c DO 4011 WHILE (T.GE.DBLE(TF)) 4009 CONTINUE IF (T.GE.DBLE(TF)) THEN CKA=CK(1)*IUPS(TE(1),T) CLA=CL(1)*IUPS(TE(1),T) DO 4012 I=1,(N-1) CKA=CKA & +((CK(I)*(TE(I+1)-T)+CK(I+1)*(T-TE(I)))/(TE(I+1)-TE(I))) & *IUPS(TE(I+1),T)*IUP(T,TE(I)) CLA=CLA & +((CL(I)*(TE(I+1)-T)+CL(I+1)*(T-TE(I)))/(TE(I+1)-TE(I))) & *IUPS(TE(I+1),T)*IUP(T,TE(I)) 4012 CONTINUE CKA=CKA+CK(N)*IUP(T,TE(N)) CLA=CLA+CL(N)*IUP(T,TE(N)) * WRITE(1,'(f8.2,2(1x,f8.4))') T,CKA,CLA jg =nmcha segini mlreel lect(icoup) = mlreel prog(1) = T prog(2) = CKA prog(3) = CLA segdes mlreel T=T-DBLE(dtemp) c4011 CONTINUE * GOTO 4009 END IF enddo segdes mlenti c WRITE(*,*)'CHAUFFAGE ENREGISTRE', t,tf goto 1000 900 continue segdes mnuag1 * la composante %m1:4 du nuage %i1 est absente ou de type incorrect call erreur(921) RETURN 1000 continue do ivar = 1,nvar1 nuavf1 = indcer(ivar) segdes nuavf1 enddo return END