Télécharger chmslv.eso

Retour à la liste

Numérotation des lignes :

  1. C CHMSLV SOURCE CHAT 05/01/12 22:00:06 5004
  2. SUBROUTINE CHMSLV(IDSCHI,SP2,ITMAX,EPS,ICALCLOG,IEM)
  3. C========================================================================
  4. C
  5. C SP ISSU DE TRIOEF ( TRSOLV)
  6. C
  7. C Modif PhM
  8. C a) Prise en charge du calcul en log de concentration
  9. c (argument ICALCLOG)
  10. C
  11. C========================================================================
  12. IMPLICIT INTEGER(I-N)
  13. IMPLICIT REAL*8(A-H,O-Z)
  14. SEGMENT IDSCHI
  15. REAL*8 GK(NYDIM),AA(NYDIM,NXDIM),FF(NZDIM,NPDIM)
  16. INTEGER IDX(NXDIM),IDY(NYDIM),IDZ(NZDIM),IDP(NPDIM),NN(6)
  17. INTEGER IDECY(NYDIM),IONZ(NXDIM)
  18. CHARACTER*32 NAME(NXDIM),NAMESP(NYDIM)
  19. ENDSEGMENT
  20. SEGMENT SP2
  21. REAL*8 GX(NXDIM),XX(NXDIM),GS(NZDIM),SS(NZDIM)
  22. REAL*8 TOT(NXDIM),TOTAQ(NXDIM),TOTFIX(NXDIM),GKS(NZDIM)
  23. REAL*8 YY(NXDIM),ZZ(NXDIM,NXDIM),CC(NYDIM),GC(NYDIM)
  24. ENDSEGMENT
  25. NXDIM=IDX(/1)
  26. NYDIM=IDY(/1)
  27. NZDIM=IDZ(/1)
  28. NPDIM=IDP(/1)
  29. ITER=0
  30.  
  31.  
  32. c write(6,*) 'ICALCLOG',ICALCLOG
  33.  
  34.  
  35. C
  36. *************************************************************************
  37. C PRISE EN COMPTE DES SOLUTIONS SOLIDES : METHODE 1
  38. N4S=0
  39. * IF(NZDIM.NE.0)THEN
  40. * I3S=NN(1)+NN(2)+NN(3)+1
  41. * I4S=NN(1)+NN(2)+NN(3)+NN(4)
  42. * DO 13 I7S=I3S,I4S
  43. * IDY7=IDY(I7S)
  44. * CALL CHIADY(IDZ,NZDIM,IDY7,ID7)
  45. * IF(ID7.NE.0)THEN
  46. * N4S=N4S+1
  47. * CALL CHMREX(IDSCHI,LGKMOD,LGKTMP,IDY7,4,5)
  48. * CALL CHMREX(IDSCHI,LGKMOD,LGKTMP,IDY7,5,4)
  49. * ENDIF
  50. * 13 CONTINUE
  51. * ENDIF
  52. **************************************************************************
  53.  
  54. NC=NN(1)+NN(2)
  55. NX=NXDIM-NN(3)-NN(4)+N4S
  56. C write(6,*)'chmslv nc nx gk',nc,nx
  57. C write(6,120)(idy(j),gk(j),j=1,nc)
  58. DO 1 J=1,NX
  59. YY(J)=0.D0
  60. 1 CONTINUE
  61.  
  62.  
  63.  
  64. 1000 CONTINUE
  65. 120 FORMAT(6(1X,I6,1PD15.6))
  66. C COMPLEXES
  67. DO 2 I=1,NC
  68. V=GK(I)
  69. * write(6,*)'chmslv idy',idy(i),'gk',gk(i)
  70. DO 3 J=1,NX
  71. V=V+AA(I,J)*GX(J)
  72. * write(6,*)'chmslv idx',idx(j)
  73. * write(6,*)'chmslv aa',aa(i,j),'gx',gx(j)
  74. 3 CONTINUE
  75. GC(I)=V
  76. C
  77. CC(I)=10.D0**GC(I)
  78. * write(6,*)'chmslv idy',IDY(I),'c',CC(I),'logc',GC(I)
  79. 2 CONTINUE
  80.  
  81.  
  82. C MOLE BALANCE
  83.  
  84. C ----------------- A REVOIR----------
  85. C IF(IADSORB.GT.0) CALL TRSURT(SP1,SP2,SURFACE)
  86. C --------------------------------------
  87. DO 201 J=1,NX
  88. V=-TOT(J)
  89. DO 200 I=1,NC
  90. V=V+AA(I,J)*CC(I)
  91. 200 CONTINUE
  92. YY(J)=V
  93. C
  94. 201 CONTINUE
  95.  
  96. C CALCUL Z (Jacobienne pour les iterations de newton)
  97. DO 301 I=1,NX
  98. DO 300 J=1,NX
  99. ZZ(I,J)=0.D0
  100. 300 CONTINUE
  101. 301 CONTINUE
  102. DO 410 J=1,NX
  103. DO 405 I=1,NC
  104. DO 400 K=1,NX
  105. c modif PhM
  106. IF (ICALCLOG.EQ.1) THEN
  107. ZZ(J,K)=ZZ(J,K)+AA(I,J)*AA(I,K)*CC(I)*LOG(10.)
  108. ELSE
  109. ZZ(J,K)=ZZ(J,K)+AA(I,J)*AA(I,K)*CC(I)/XX(K)
  110. END IF
  111. c modif PhM
  112. 400 CONTINUE
  113. 405 CONTINUE
  114. C write(6,*)'chmslv iter zz',iter,'IDX ',IDX(J)
  115. c write(6,120)(idx(k),ZZ(J,K),K=1,NX)
  116. 410 CONTINUE
  117. C
  118. C
  119. C ----------------- A REVOIR----------
  120. C IF (IADSORB.GT.0) CALL TRSURZ(SP1,SP2,SURFACE)
  121. C ----------------------------------------------
  122.  
  123. C CONVERGENCE TEST
  124.  
  125. DO 800 J=1,NX
  126. TJ = TOT(J)
  127. VMAX=ABS(TJ)
  128. DO 810 I=1,NC
  129. IF(ABS(AA(I,J)*CC(I)).GE.VMAX) THEN
  130. VMAX=ABS(AA(I,J)*CC(I))
  131. ENDIF
  132. 810 CONTINUE
  133. YJ = YY(J)
  134. C pour test et impression
  135. C idximp=idx(j)
  136. C WRITE(6,*)' CHMSLV J YJ VMAX EPS ',J,YJ,VMAX,EPS
  137. c write(6,*) ' J YJ/VMAX ',J,ABS(YJ)/VMAX
  138. IF (ABS(YJ)/VMAX.GT.EPS) GOTO 840
  139. 800 CONTINUE
  140.  
  141.  
  142. c write(6,*) 'convergence en ', iter, ' etapes'
  143. c write(6,120) (idx(j),gx(j),j=1,nx)
  144.  
  145. RETURN
  146.  
  147. 840 CONTINUE
  148.  
  149. ITER=ITER+1
  150. IF(ITER.GT.ITMAX) THEN
  151. C WRITE(6,*)' CHMSLV YJ VMAX IDX TJ',YJ,VMAX,IDXIMP,TJ
  152. C WRITE(6,*) ' PROBLEME DE CONVERGENCE EN CHIMIE AU NOEUD ',NI
  153. C WRITE(6,*) ' ITER: ',ITER,'> ITMAX: ',ITMAX
  154. C WRITE(6,*) ' L ETAT DU PROBLEME EST LE SUIVANT: '
  155. C WRITE(6,*) 'EXECUTION CHIMIE TERMINEE ** ERREUR 7'
  156. IEM = 1
  157. RETURN
  158. ENDIF
  159.  
  160. C WRITE(6,*) 'CONVERGENCE TEST OK'
  161. C
  162. C ITERATE
  163.  
  164. C Resolution du systeme lineaire
  165. CALL CHMSMQ(ZZ,YY,NX,NXDIM,IEM)
  166. C write(6,*)'chmslv xx',iter
  167. c write(6,120)(idx(j),yy(j),j=1,nx)
  168.  
  169.  
  170.  
  171. IF(IEM.NE.0)RETURN
  172. C
  173.  
  174. C
  175. C Calcul de la concentration a l'etape l
  176. C
  177. DO 500 J=1,NX
  178. c modif PhM
  179. c Deux moyen de calculer la concentration :
  180. c C^l = 10^{LC^{l-1} - dLC^{l}} ou
  181. c C^l = C^{l-1} * 10^{-dLC^{l}}
  182. c XX(J) = 10.D0**(GX(J))
  183. IF (ICALCLOG.EQ.1) THEN
  184. GX(J) = GX(J)-YY(J)
  185. XX(J) = XX(J) * 10.D0**(-YY(J))
  186. ELSE
  187. XX(J)=XX(J)-YY(J)
  188. IF (XX(J).LE.0.D0) THEN
  189. XX(J)=(XX(J)+YY(J))/10.D0
  190. ENDIF
  191. c modif PhM
  192. c IF (ABS(XX(J)).LT.1.D-60) XX(J)=1.D-60
  193. IF (ABS(XX(J)).LT.1.D-100) XX(J)=1.D-100
  194. c modif PhM
  195. C WRITE(6,510) XX(J)
  196. GX(J)=LOG10(XX(J))
  197. END IF
  198. C modif PhM
  199.  
  200. * write(6,*)'chmslv idx',idx(j),'gx',gx(j)
  201.  
  202. 500 CONTINUE
  203. 510 FORMAT ('XX(J): ', 1PE19.10)
  204. GO TO 1000
  205. END
  206.  
  207.  
  208.  
  209.  
  210.  
  211.  
  212.  
  213.  
  214.  
  215.  
  216.  
  217.  
  218.  
  219.  

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