C CHMSLV SOURCE CHAT 05/01/12 22:00:06 5004 SUBROUTINE CHMSLV(IDSCHI,SP2,ITMAX,EPS,ICALCLOG,IEM) C======================================================================== C C SP ISSU DE TRIOEF ( TRSOLV) C C Modif PhM C a) Prise en charge du calcul en log de concentration c (argument ICALCLOG) C C======================================================================== IMPLICIT INTEGER(I-N) IMPLICIT REAL*8(A-H,O-Z) SEGMENT IDSCHI REAL*8 GK(NYDIM),AA(NYDIM,NXDIM),FF(NZDIM,NPDIM) INTEGER IDX(NXDIM),IDY(NYDIM),IDZ(NZDIM),IDP(NPDIM),NN(6) INTEGER IDECY(NYDIM),IONZ(NXDIM) CHARACTER*32 NAME(NXDIM),NAMESP(NYDIM) ENDSEGMENT SEGMENT SP2 REAL*8 GX(NXDIM),XX(NXDIM),GS(NZDIM),SS(NZDIM) REAL*8 TOT(NXDIM),TOTAQ(NXDIM),TOTFIX(NXDIM),GKS(NZDIM) REAL*8 YY(NXDIM),ZZ(NXDIM,NXDIM),CC(NYDIM),GC(NYDIM) ENDSEGMENT NXDIM=IDX(/1) NYDIM=IDY(/1) NZDIM=IDZ(/1) NPDIM=IDP(/1) ITER=0 c write(6,*) 'ICALCLOG',ICALCLOG C ************************************************************************* C PRISE EN COMPTE DES SOLUTIONS SOLIDES : METHODE 1 N4S=0 * IF(NZDIM.NE.0)THEN * I3S=NN(1)+NN(2)+NN(3)+1 * I4S=NN(1)+NN(2)+NN(3)+NN(4) * DO 13 I7S=I3S,I4S * IDY7=IDY(I7S) * CALL CHIADY(IDZ,NZDIM,IDY7,ID7) * IF(ID7.NE.0)THEN * N4S=N4S+1 * CALL CHMREX(IDSCHI,LGKMOD,LGKTMP,IDY7,4,5) * CALL CHMREX(IDSCHI,LGKMOD,LGKTMP,IDY7,5,4) * ENDIF * 13 CONTINUE * ENDIF ************************************************************************** NC=NN(1)+NN(2) NX=NXDIM-NN(3)-NN(4)+N4S C write(6,*)'chmslv nc nx gk',nc,nx C write(6,120)(idy(j),gk(j),j=1,nc) DO 1 J=1,NX YY(J)=0.D0 1 CONTINUE 1000 CONTINUE 120 FORMAT(6(1X,I6,1PD15.6)) C COMPLEXES DO 2 I=1,NC V=GK(I) * write(6,*)'chmslv idy',idy(i),'gk',gk(i) DO 3 J=1,NX V=V+AA(I,J)*GX(J) * write(6,*)'chmslv idx',idx(j) * write(6,*)'chmslv aa',aa(i,j),'gx',gx(j) 3 CONTINUE GC(I)=V C CC(I)=10.D0**GC(I) * write(6,*)'chmslv idy',IDY(I),'c',CC(I),'logc',GC(I) 2 CONTINUE C MOLE BALANCE C ----------------- A REVOIR---------- C IF(IADSORB.GT.0) CALL TRSURT(SP1,SP2,SURFACE) C -------------------------------------- DO 201 J=1,NX V=-TOT(J) DO 200 I=1,NC V=V+AA(I,J)*CC(I) 200 CONTINUE YY(J)=V C 201 CONTINUE C CALCUL Z (Jacobienne pour les iterations de newton) DO 301 I=1,NX DO 300 J=1,NX ZZ(I,J)=0.D0 300 CONTINUE 301 CONTINUE DO 410 J=1,NX DO 405 I=1,NC DO 400 K=1,NX c modif PhM IF (ICALCLOG.EQ.1) THEN ZZ(J,K)=ZZ(J,K)+AA(I,J)*AA(I,K)*CC(I)*LOG(10.) ELSE ZZ(J,K)=ZZ(J,K)+AA(I,J)*AA(I,K)*CC(I)/XX(K) END IF c modif PhM 400 CONTINUE 405 CONTINUE C write(6,*)'chmslv iter zz',iter,'IDX ',IDX(J) c write(6,120)(idx(k),ZZ(J,K),K=1,NX) 410 CONTINUE C C C ----------------- A REVOIR---------- C IF (IADSORB.GT.0) CALL TRSURZ(SP1,SP2,SURFACE) C ---------------------------------------------- C CONVERGENCE TEST DO 800 J=1,NX TJ = TOT(J) VMAX=ABS(TJ) DO 810 I=1,NC IF(ABS(AA(I,J)*CC(I)).GE.VMAX) THEN VMAX=ABS(AA(I,J)*CC(I)) ENDIF 810 CONTINUE YJ = YY(J) C pour test et impression C idximp=idx(j) C WRITE(6,*)' CHMSLV J YJ VMAX EPS ',J,YJ,VMAX,EPS c write(6,*) ' J YJ/VMAX ',J,ABS(YJ)/VMAX IF (ABS(YJ)/VMAX.GT.EPS) GOTO 840 800 CONTINUE c write(6,*) 'convergence en ', iter, ' etapes' c write(6,120) (idx(j),gx(j),j=1,nx) RETURN 840 CONTINUE ITER=ITER+1 IF(ITER.GT.ITMAX) THEN C WRITE(6,*)' CHMSLV YJ VMAX IDX TJ',YJ,VMAX,IDXIMP,TJ C WRITE(6,*) ' PROBLEME DE CONVERGENCE EN CHIMIE AU NOEUD ',NI C WRITE(6,*) ' ITER: ',ITER,'> ITMAX: ',ITMAX C WRITE(6,*) ' L ETAT DU PROBLEME EST LE SUIVANT: ' C WRITE(6,*) 'EXECUTION CHIMIE TERMINEE ** ERREUR 7' IEM = 1 RETURN ENDIF C WRITE(6,*) 'CONVERGENCE TEST OK' C C ITERATE C Resolution du systeme lineaire CALL CHMSMQ(ZZ,YY,NX,NXDIM,IEM) C write(6,*)'chmslv xx',iter c write(6,120)(idx(j),yy(j),j=1,nx) IF(IEM.NE.0)RETURN C C C Calcul de la concentration a l'etape l C DO 500 J=1,NX c modif PhM c Deux moyen de calculer la concentration : c C^l = 10^{LC^{l-1} - dLC^{l}} ou c C^l = C^{l-1} * 10^{-dLC^{l}} c XX(J) = 10.D0**(GX(J)) IF (ICALCLOG.EQ.1) THEN GX(J) = GX(J)-YY(J) XX(J) = XX(J) * 10.D0**(-YY(J)) ELSE XX(J)=XX(J)-YY(J) IF (XX(J).LE.0.D0) THEN XX(J)=(XX(J)+YY(J))/10.D0 ENDIF c modif PhM c IF (ABS(XX(J)).LT.1.D-60) XX(J)=1.D-60 IF (ABS(XX(J)).LT.1.D-100) XX(J)=1.D-100 c modif PhM C WRITE(6,510) XX(J) GX(J)=LOG10(XX(J)) END IF C modif PhM * write(6,*)'chmslv idx',idx(j),'gx',gx(j) 500 CONTINUE 510 FORMAT ('XX(J): ', 1PE19.10) GO TO 1000 END