1 : C SSCH SOURCE CHAT 06/03/29 21:34:55 5360 2 : subroutine ssch 3 : 4 : ccccccccccccccccccccccccccccccccccccccccccc 5 : c calcul des terme source si transport par espece 6 : c 7 : c Ri = ssch tabchi1 Ai SOLU dsol 8 : c 9 : c Ri terme source nesp chpoi 10 : c Ai variation sur les composantes fixees ncomp chpoi 11 : c SOLU concentrations des especes nesp chpoi 12 : c dsol c1-c2 nesp chpoi 13 : c 14 : cccccccccccccccccccccccccccc 15 : IMPLICIT INTEGER(I-N) 16 : IMPLICIT REAL*8(A-H,O-Z) 17 : -inc CCOPTIO 18 : -inc SMLENTI 19 : -inc SMLREEL 20 : -inc SMTABLE 21 : -inc SMLMOTS 22 : -inc SMCHPOI 23 : -inc SMELEME 24 : -inc SMCOORD 25 : 26 : pointeur mlaa.mlreel,mlogk.mlreel,mlff.mlreel 27 : pointeur mlidx.mlenti,mlidy.mlenti,mlnn.mlenti,mldecy.mlenti 28 : pointeur mlionz.mlenti,mlsolu.mlenti 29 : pointeur mlidp.mlenti,mlidz.mlenti 30 : pointeur mlname.mlmots,mlnesp.mlmots 31 : pointeur izt.mchpoi,izs.mchpoi,mchrij.mchpoi,izd.mchpoi 32 : pointeur ipt.mpoval,ips.mpoval,ichrij.mpoval,ipd.mpoval 33 : 34 : segment idschi 35 : real*8 gk(nydim),aa(nydim,nxdim),ff(nzdim,npdim) 36 : integer idx(nxdim),idy(nydim),idz(nzdim),idp(npdim),nn(6) 37 : integer idecy(nydim),ionz(nxdim) 38 : character*32 name(nxdim),namesp(nydim) 39 : endsegment 40 : segment sp2 41 : real*8 gx(nxdim),xx(nxdim),gs(nzdim),ss(nzdim) 42 : real*8 tot(nxdim),totaq(nxdim),totfix(nxdim),gks(nzdim) 43 : real*8 yy(nxdim),zz(nxdim,nxdim),cc(nydim),gc(nydim) 44 : endsegment 45 : segment iptidx 46 : integer itdx(nxdim) 47 : endsegment 48 : pointeur idxtot.iptidx,idyc.iptidx 49 : segment iztra 50 : integer itab(nval),ktab(nn12) 51 : endsegment 52 : 53 : character*8 ltyp 54 : character*4 nomtot 55 : 56 : 57 : 58 : 59 : 60 : 61 : c ** lecture table issue de chi1 62 : c 62 : call chmdeb(mlaa,mlogk,mlff,mlidx,mlidy,mlidz,mlidp,mlnn,mldecy, 63 : * mlname,mlionz,itiden,itredo,itempe,mlnesp) 64 : 65 : if(ierr.ne.0) return 66 : 67 : c ** recuperation des chpoints 68 : c 69 : 70 : ltyp = 'CHPOINT ' 71 : iretou = 0 72 : icod = 1 72 : call lirobj( ltyp , iret, icod,iretou) 73 : if(iretou.eq.0)return 74 : izt = iret 75 : segact izt 76 : 77 : iretou = 0 78 : icod = 1 78 : call lirobj( ltyp , iret, icod,iretou) 79 : if(iretou.eq.0)return 80 : izs = iret 81 : segact izs 82 : 83 : iretou = 0 84 : icod = 1 84 : call lirobj( ltyp , iret, icod,iretou) 85 : if(iretou.eq.0)return 86 : izd = iret 87 : segact izd 88 : 89 : 90 : c ** echange eventuel A1 SOLU 91 : c 92 : 93 : c ** coherence des support de Ai et SOLU 94 : c 95 : 96 : nsoupo = izt.ipchp(/1) 97 : if( nsoupo.ne.1) then 98 : moterr(1:8)= 'CHAMPOIN' 98 : call erreur(132) 99 : segdes izt 100 : endif 101 : msoupo = izt.ipchp(1) 102 : segact msoupo 103 : meleme=igeoc 104 : indic=1 105 : nbcomp =-1 106 : nomtot=' ' 106 : call quepoi(izs,meleme,indiq,nbcomp,nomtot) 107 : if(indiq.lt.0) then 107 : call erreur(22) 108 : endif 109 : if(ierr.ne.0)return 110 : 110 : call quepoi(izd,meleme,indiq,nbcomp,nomtot) 111 : if(indiq.lt.0) then 111 : call erreur(22) 112 : endif 113 : if(ierr.ne.0)return 114 : 115 : c ** lecture de la table iden les segments reviennent actifs 116 : c 116 : call chmide(itiden,mlcomp,mlsolu,mmsolu,mlprec,mmprec,mlsurf, 117 : * mmsurf,mltyp3,mmtyp3,mltyp6,mmtyp6,mlparf,mlreac,mlimmo, 118 : * mlpole,mmpole,mlsoso,mmsoso,limp3) 119 : if(ierr.ne.0) return 120 : 121 : segact mlaa,mlogk,mlidx,mlidy,mlnn,mldecy,mlname,mlionz,mlnesp 122 : segact mlff,mlidz,mlidp 123 : nxdim=mlidx.lect(/1) 124 : nydim=mlidy.lect(/1) 125 : nzdim=mlidz.lect(/1) 126 : npdim=mlidp.lect(/1) 127 : segini idschi 128 : segini sp2 129 : 129 : call chmids(mlaa,mlogk,mlff,mlidx,mlidy,mlidz,mlidp,mlnn,mldecy, 130 : * mlname,mlionz,idschi,mlnesp) 131 : 132 : nn12 = mlsolu.lect(/1) 133 : n1n2 = nn(1) +nn(2) 134 : 135 : nval = nn12 - nxdim 136 : segini iztra 137 : 138 : 139 : moterr(5:8)='TOT ' 139 : call chmlst (izt,mlidx,idxtot,ipt ) 140 : 141 : moterr(5:8)='SOLU' 141 : call chmmst (izs,mmsolu,idyc ,ips ) 142 : 143 : segact izs 144 : msoupo = izs.ipchp(1) 145 : segact msoupo 146 : ips = ipoval 147 : segact ips 148 : segdes msoupo 149 : 150 : segact izd 151 : msoupo = izd.ipchp(1) 152 : segact msoupo 153 : ipd = ipoval 154 : segact ipd 155 : segdes msoupo 156 : 157 : npn = ipt.vpocha(/1) 158 : 159 : 160 : c ** creation chpoi de resultats 161 : c 161 : call chmcrc(mmsolu,meleme,npn,mchrij,ichrij) 162 : * write(6,*)' chmcrc ' 163 : 164 : do 100 ii =1,npn 165 : 166 : c ** chargement de idschi 167 : c 167 : call chmids(mlaa,mlogk,mlff,mlidx,mlidy,mlidz,mlidp,mlnn,mldecy, 168 : * mlname,mlionz,idschi,mlnesp) 169 : 170 : do j=1,nxdim 171 : ipg = idxtot.itdx(j) 172 : tot(j)=ipt.vpocha(ii,ipg ) 173 : enddo 174 : 175 : ites = 0 176 : do j=1,nxdim 177 : if(abs(tot(j)).gt.1.d-31) then 178 : ites = ites + 1 179 : endif 180 : enddo 181 : 182 : if(ites.gt.0) then 183 : do j=1, nn12 184 : jpg = idyc.itdx(j) 185 : cc(j)=ips.vpocha(ii, jpg ) 186 : enddo 187 : 188 : do j =1,nxdim 189 : idxt = idx(j) 190 : do k =1,nydim 191 : if(idy(k).eq.idxt) then 192 : xx(j) = cc(k) 193 : go to 10 194 : endif 195 : enddo 196 : 10 continue 197 : enddo 198 : 199 : do i=1,nxdim 200 : do j=1,nxdim 201 : zz(i,j)=0.d0 202 : do l=1,nn12 203 : idyt = mlsolu.lect(l) 204 : do ll =1,n1n2 205 : if(idy(ll).eq.idyt) then 206 : zz(i,j)= zz(i,j) + aa(ll,i) * aa(ll,j) * cc( l) /xx(j) 207 : go to 20 208 : endif 209 : enddo 210 : 20 continue 211 : enddo 212 : enddo 213 : enddo 214 : 215 : 216 : 217 : iem = 0 217 : call chmsmq(zz,tot,nxdim,nxdim,iem) 218 : 219 : do l = 1,nn12 220 : gc(l) = 0.d0 221 : idyt = mlsolu.lect(l) 222 : do ll =1,n1n2 223 : if(idy(ll).eq.idyt) then 224 : do j =1,nxdim 225 : gc(l)=gc(l) + aa(ll,j) * cc(l) / xx(j) * tot(j) 226 : enddo 227 : go to 40 228 : endif 229 : enddo 230 : 40 continue 231 : enddo 232 : 233 : kk = 0 234 : do k=1,nn12 235 : ichrij.vpocha(ii, k)= gc(k) 236 : 30 continue 237 : enddo 238 : 239 : else 240 : * do k = 1,nn12 241 : * ichrij . vpocha(ii,k) = ipd . vpocha(ii,k) 242 : * enddo 243 : * endif 244 : 245 : * ioldf = 1 246 : * if(ioldf.eq.156) then 247 : 248 : 249 : if(nval.ge.nxdim) then 250 : do i = 1 ,nval 251 : rr = ipd.vpocha(ii,i) 252 : idyt = mlsolu.lect(i) 253 : do il = 1,n1n2 254 : if(idy(il).eq.idyt)then 255 : do j=1,nxdim 256 : tot(j) = tot(j) - aa(il,j) * rr 257 : enddo 258 : go to 50 259 : endif 260 : enddo 261 : 50 continue 262 : enddo 263 : 264 : do i = nval + 1 , nn12 265 : idyt = mlsolu.lect(i) 266 : ik = i - nval 267 : do il = 1,n1n2 268 : if(idy(il).eq.idyt)then 269 : do j=1,nxdim 270 : zz(j ,ik) = aa(il,j) 271 : enddo 272 : go to 60 273 : endif 274 : enddo 275 : 60 continue 276 : enddo 277 : 277 : call chmsmq(zz,tot,nxdim,nxdim,iem) 278 : 279 : do k =1,nval 280 : ichrij.vpocha(ii,k) = ipd.vpocha(ii,k) 281 : enddo 282 : do k = nval+1,nn12 283 : ik = k - nval 284 : ichrij.vpocha(ii,k) = tot ( ik) 285 : enddo 286 : else 287 : 288 : do i=1,nval 289 : itab(i)=0 290 : enddo 291 : do i=1,nn12 292 : ktab(i)=0 293 : enddo 294 : do j =1,nxdim 295 : gx(j) = ips.vpocha(ii,j) 296 : gc(j) = ipd.vpocha(ii,j) 297 : enddo 298 : do j=1,nval 299 : vmi = 100. 300 : do k =1,nxdim 301 : if(abs(gx(k)).lt.vmi.and.abs(gx(k)).gt.1.d-37) then 302 : vmi = gx(k) 303 : im = k 304 : endif 305 : enddo 306 : gx(im) = 100. 307 : itab(j) = im 308 : idyt = mlsolu.lect(im) 309 : rr = ipd.vpocha(ii,im) 310 : do il = 1,n1n2 311 : if(idy(il).eq.idyt)then 312 : do kj=1,nxdim 313 : tot(kj) = tot(kj) - aa(il,kj) * rr 314 : enddo 315 : go to 70 316 : endif 317 : enddo 318 : 70 continue 319 : enddo 320 : do ik=1,nval 321 : ktab(itab(ik))=ik 322 : enddo 323 : 324 : iki = 0 325 : do i = 1 , nn12 326 : if(ktab(i).eq.0) then 327 : idyt = mlsolu.lect(i) 328 : iki = iki + 1 329 : do il = 1,n1n2 330 : if(idy(il).eq.idyt)then 331 : do j=1,nxdim 332 : zz(j,iki) = aa(il,j) 333 : enddo 334 : go to 80 335 : endif 336 : enddo 337 : 80 continue 338 : endif 339 : enddo 340 : 340 : call chmsmq(zz,tot,nxdim,nxdim,iem) 341 : 342 : kk = 0 343 : do k =1,nn12 344 : if(ktab(k).ne.0) then 345 : ki = itab(ktab(k)) 346 : ichrij.vpocha(ii,k) = ipd.vpocha(ii,ki) 347 : else 348 : kk = kk+ 1 349 : 350 : 351 : ichrij.vpocha(ii,k) = tot ( kk) 352 : endif 353 : enddo 354 : 355 : endif 356 : 357 : 358 : endif 359 : 360 : 100 continue 361 : 362 : c ** menage 363 : c 363 : call ecrobj('CHPOINT ',mchrij) 364 : mchpoi = mchrij 365 : msoupo = ipchp(1) 366 : mpoval = ipoval 367 : segdes mchpoi,msoupo,mpoval 368 : 369 : 370 : segsup idschi 371 : segsup sp2 372 : segsup idxtot 373 : segsup iztra 374 : 375 : segdes mlaa,mlogk,mlidx,mlidy,mlnn,mldecy,mlname,mlionz,mlnesp 376 : 377 : segdes mlsolu 378 : segdes meleme 379 : 380 : segact izt 381 : msoupo = izt.ipchp(1) 382 : segdes msoupo 383 : segact izs 384 : msoupo = izs.ipchp(1) 385 : segdes msoupo 386 : segact izd 387 : msoupo = izd.ipchp(1) 388 : segdes msoupo 389 : segdes ipt,ips,ipd 390 : segdes izt,izs,izd 391 : 392 : if(mlimmo.ne.0) then 393 : mlenti=mlimmo 394 : segdes mlenti 395 : endif 396 : if(mlreac.ne.0) then 397 : mlenti=mlreac 398 : segdes mlenti 399 : endif 400 : if(mlparf.ne.0) then 401 : mlenti= mlparf 402 : segdes mlenti 403 : endif 404 : if(mltyp6.ne.0) then 405 : mlenti=mltyp6 406 : segdes mlenti 407 : endif 408 : if(mltyp3.ne.0) then 409 : mlenti=mltyp3 410 : segdes mlenti 411 : endif 412 : if(mlsurf.ne.0) then 413 : mlenti=mlsurf 414 : segdes mlenti 415 : endif 416 : if(mlprec.ne.0) then 417 : mlenti=mlprec 418 : segdes mlenti 419 : endif 420 : if(mlcomp.ne.0) then 421 : mlenti=mlcomp 422 : segdes mlenti 423 : endif 424 : if(mmtyp6.ne.0) then 425 : mlmots=mmtyp6 426 : segdes mlmots 427 : endif 428 : if(mmtyp3.ne.0) then 429 : mlmots=mmtyp3 430 : segdes mlmots 431 : endif 432 : if(mmsurf.ne.0) then 433 : mlmots=mmsurf 434 : segdes mlmots 435 : endif 436 : if(mmprec.ne.0) then 437 : mlmots=mmprec 438 : segdes mlmots 439 : endif 440 : 441 : mtab2 =itiden 442 : segdes mtab2 443 : 444 : 445 : 446 : return 447 : end 448 : 449 : 450 : 451 :
© Cast3M 2003 - All rights reserved.
Disclaimer