Télécharger la plaquette Cast3M


fr en

Télécharger ssch.eso

Retour à la liste

Numérotation des lignes : 

ssch

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

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