!------------------------------------------------------------------------------- ! This file is part of the Code_Saturne Kernel, element of the ! Code_Saturne CFD tool. ! Copyright (C) 1998-2009 EDF S.A., France ! contact: saturne-support@edf.fr ! The Code_Saturne Kernel is free software; you can redistribute it ! and/or modify it under the terms of the GNU General Public License ! as published by the Free Software Foundation; either version 2 of ! the License, or (at your option) any later version. ! The Code_Saturne Kernel is distributed in the hope that it will be ! useful, but WITHOUT ANY WARRANTY; without even the implied warranty ! of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the ! GNU General Public License for more details. ! You should have received a copy of the GNU General Public License ! along with the Code_Saturne Kernel; if not, write to the ! Free Software Foundation, Inc., ! 51 Franklin St, Fifth Floor, ! Boston, MA 02110-1301 USA !------------------------------------------------------------------------------- subroutine raypun & !================ ( idbia0 , idbra0 , & ndim , ncelet , ncel , nfac , nfabor , nfml , nprfml , & nnod , lndfac , lndfbr , ncelbr , & nvar , nscal , nphas , iphas , & nideve , nrdeve , nituse , nrtuse , & ifacel , ifabor , ifmfbr , ifmcel , iprfml , itypfb , & ipnfac , nodfac , ipnfbr , nodfbr , & idevel , ituser , ia , & xyzcen , surfac , surfbo , cdgfac , cdgfbo , xyznod , volume , & dt , rtp , rtpa , propce , propfa , propfb , & coefa , coefb , & cofrua , cofrub , & flurds , flurdb , & dtr , viscf , viscb , & dam , xam , & drtp , smbrs , rovsdt , & theta4 , thetaa , sa , & qx , qy , qz , & qincid , eps , tparoi , & w1 , w2 , w3 , w4 , w5 , & w6 , w7 , w8 , w9 , ckmel , & rdevel , rtuser , ra ) !=============================================================================== ! FONCTION : ! ---------- ! SOUS-PROGRAMME DU MODULE DE RAYONNEMENT : ! ----------------------------------------- ! CALCUL DES FLUX ET DU TERME SOURCE RADIATIFS ! AVEC L'APPROXIMATION P-1 !------------------------------------------------------------------------------- !ARGU ARGUMENTS !__________________.____._____.________________________________________________. ! name !type!mode ! role ! !__________________!____!_____!________________________________________________! ! idbia0 ! i ! <-- ! number of first free position in ia ! ! idbra0 ! i ! <-- ! number of first free position in ra ! ! ndim ! i ! <-- ! spatial dimension ! ! ncelet ! i ! <-- ! number of extended (real + ghost) cells ! ! ncel ! i ! <-- ! number of cells ! ! nfac ! i ! <-- ! number of interior faces ! ! nfabor ! i ! <-- ! number of boundary faces ! ! nfml ! i ! <-- ! number of families (group classes) ! ! nprfml ! i ! <-- ! number of properties per family (group class) ! ! nnod ! i ! <-- ! number of vertices ! ! lndfac ! i ! <-- ! size of nodfac indexed array ! ! lndfbr ! i ! <-- ! size of nodfbr indexed array ! ! ncelbr ! i ! <-- ! number of cells with faces on boundary ! ! nvar ! i ! <-- ! total number of variables ! ! nscal ! i ! <-- ! total number of scalars ! ! nphas ! i ! <-- ! number of phases ! ! iphas ! i ! --> ! phase number ! ! nideve, nrdeve ! i ! <-- ! sizes of idevel and rdevel arrays ! ! nituse, nrtuse ! i ! <-- ! sizes of ituser and rtuser arrays ! ! ifacel(2, nfac) ! ia ! <-- ! interior faces -> cells connectivity ! ! ifabor(nfabor) ! ia ! <-- ! boundary faces -> cells connectivity ! ! ifmfbr(nfabor) ! ia ! <-- ! boundary face family numbers ! ! ifmcel(ncelet) ! ia ! <-- ! cell family numbers ! ! iprfml ! ia ! <-- ! property numbers per family ! ! (nfml, nprfml) ! ! ! ! ! itypfb ! ia ! <-- ! boundary face types ! ! (nfabor, nphas) ! ! ! ! ! ipnfac(nfac+1) ! ia ! <-- ! interior faces -> vertices index (optional) ! ! nodfac(lndfac) ! ia ! <-- ! interior faces -> vertices list (optional) ! ! ipnfbr(nfabor+1) ! ia ! <-- ! boundary faces -> vertices index (optional) ! ! nodfbr(lndfbr) ! ia ! <-- ! boundary faces -> vertices list (optional) ! ! idevel(nideve) ! ia ! <-> ! integer work array for temporary development ! ! ituser(nituse) ! ia ! <-> ! user-reserved integer work array ! ! ia(*) ! ia ! --- ! main integer work array ! ! xyzcen ! ra ! <-- ! cell centers ! ! (ndim, ncelet) ! ! ! ! ! surfac ! ra ! <-- ! interior faces surface vectors ! ! (ndim, nfac) ! ! ! ! ! surfbo ! ra ! <-- ! boundary faces surface vectors ! ! (ndim, nfabor) ! ! ! ! ! cdgfac ! ra ! <-- ! interior faces centers of gravity ! ! (ndim, nfac) ! ! ! ! ! cdgfbo ! ra ! <-- ! boundary faces centers of gravity ! ! (ndim, nfabor) ! ! ! ! ! xyznod ! ra ! <-- ! vertex coordinates (optional) ! ! (ndim, nnod) ! ! ! ! ! volume(ncelet) ! ra ! <-- ! cell volumes ! ! dt(ncelet) ! ra ! <-- ! time step (per cell) ! ! rtp, rtpa ! ra ! <-- ! calculated variables at cell centers ! ! (ncelet, *) ! ! ! (at current and previous time steps) ! ! propce(ncelet, *)! ra ! <-- ! physical properties at cell centers ! ! propfa(nfac, *) ! ra ! <-- ! physical properties at interior face centers ! ! propfb(nfabor, *)! ra ! <-- ! physical properties at boundary face centers ! ! coefa, coefb ! ra ! <-- ! boundary conditions ! ! (nfabor, *) ! ! ! ! ! cofrua,cofrub ! tr ! --- ! conditions aux limites aux ! !(nfabor) ! ! ! faces de bord pour la luminance ! ! flurds,flurdb ! tr ! --- ! pseudo flux de masse (faces internes ! !(nfac)(nfabor) ! ! ! et faces de bord ) ! ! dtr(ncelet) ! tr ! --- ! dt*cdtvar ! ! viscf(nfac) ! tr ! --- ! visc*surface/dist aux faces internes ! ! viscb(nfabor ! tr ! --- ! visc*surface/dist aux faces de bord ! ! dam(ncelet ! tr ! --- ! tableau de travail pour matrice ! ! xam(nfac,*) ! tr ! --- ! tableau de travail pour matrice ! ! drtp(ncelet ! tr ! --- ! tableau de travail pour increment ! ! smbrs(ncelet ! tr ! --- ! tableau de travail pour sec mem ! ! rovsdt(ncelet ! tr ! --- ! tableau de travail pour terme instat ! ! theta4(ncelet ! tr ! --- ! pseudo temperature radiative ! ! thetaa(ncelet ! tr ! --- ! pseudo temp rar pdt precedent (nulle) ! ! sa (ncelet) ! tr ! --> ! part d'absorption du terme source rad ! ! qxqyqz(ncelet ! tr ! --> ! composante du vecteur densite de flux ! ! ! ! ! radiatif explicite ! ! qincid(nfabor ! tr ! --> ! densite de flux radiatif aux bords ! ! eps (nfabor) ! tr ! <-- ! emissivite des facettes de bord ! ! tparoi(nfabor ! tr ! <-- ! temperature de paroi en kelvin ! ! w1...9(ncelet ! tr ! --- ! tableau de travail ! ! ckmel(ncelet) ! tr ! <-- ! coeff d'absorption du melange ! ! ! ! ! gaz-particules de charbon ! ! rdevel(nrdeve) ! ra ! <-> ! real work array for temporary development ! ! rtuser(nrtuse) ! ra ! <-> ! user-reserved real work array ! ! ra(*) ! ra ! --- ! main real work array ! !__________________!____!_____!________________________________________________! ! TYPE : E (ENTIER), R (REEL), A (ALPHANUMERIQUE), T (TABLEAU) ! L (LOGIQUE) .. ET TYPES COMPOSES (EX : TR TABLEAU REEL) ! MODE : <-- donnee, --> resultat, <-> Donnee modifiee ! --- tableau de travail !------------------------------------------------------------------------------- !=============================================================================== implicit none !=============================================================================== ! Common blocks !=============================================================================== include "paramx.h" include "numvar.h" include "entsor.h" include "optcal.h" include "cstphy.h" include "cstnum.h" include "pointe.h" include "ppppar.h" include "ppthch.h" include "cpincl.h" include "ppincl.h" include "radiat.h" include "parall.h" include "period.h" !=============================================================================== ! Arguments integer idbia0 , idbra0 integer ndim , ncelet , ncel , nfac , nfabor integer nfml , nprfml integer nnod , lndfac , lndfbr , ncelbr integer nvar , nscal , nphas , iphas integer nideve , nrdeve , nituse , nrtuse integer ifacel(2,nfac) , ifabor(nfabor) integer ifmfbr(nfabor) , ifmcel(ncelet) integer iprfml(nfml,nprfml) , itypfb(nfabor,nphas) integer ipnfac(nfac+1), nodfac(lndfac) integer ipnfbr(nfabor+1), nodfbr(lndfbr) integer idevel(nideve), ituser(nituse) integer ia(*) double precision xyzcen(ndim,ncelet) double precision surfac(ndim,nfac), surfbo(ndim,nfabor) double precision cdgfac(ndim,nfac), cdgfbo(ndim,nfabor) double precision xyznod(ndim,nnod), volume(ncelet) double precision dt(ncelet), rtp(ncelet,*), rtpa(ncelet,*) double precision propce(ncelet,*) double precision propfa(nfac,*), propfb(nfabor,*) double precision coefa(nfabor,*), coefb(nfabor,*) double precision cofrua(nfabor), cofrub(nfabor) double precision flurds(nfac), flurdb(nfabor) double precision dtr(ncelet) double precision viscf(nfac), viscb(nfabor) double precision dam(ncelet), xam(nfac,2) double precision drtp(ncelet), smbrs(ncelet) double precision rovsdt(ncelet) double precision theta4(ncelet), thetaa(ncelet) double precision sa(ncelet) double precision qx(ncelet), qy(ncelet), qz(ncelet) double precision qincid(nfabor), tparoi(nfabor), eps(nfabor) double precision w1(ncelet), w2(ncelet), w3(ncelet) double precision w4(ncelet), w5(ncelet), w6(ncelet) double precision w7(ncelet), w8(ncelet), w9(ncelet) double precision ckmel(ncelet) double precision rdevel(nrdeve), rtuser(nrtuse), ra(*) ! Local variables character*80 cnom integer idebia, idebra integer ifac , iel , fel , numpar ! ajout de fel, numpar integer iconv1, idiff1, ndirc1, ireso1 integer nitmap, nswrsp, nswrgp, iwarnp integer imgr1 , imligp, ircflp, ischcp, isstpp, iescap integer ncymap, nitmgp integer inum integer idtva0, ivar0 integer inc, iccocg, iphydp integer idimte , itenso integer visible ! ajout de visible double precision xyzsym(3) ! ajout de xyzsym(3) double precision xpt1 , ypt1 , zpt1 ! ajout de xpt1, ypt1, zpt1 double precision xpt2 , ypt2 , zpt2 ! ajout de xpt2, ypt2, zpt2 double precision xpt3 , ypt3 , zpt3 ! ajout de xpt3, ypt3, zpt3 double precision vectux, vectuy, vectvy, vectvz ! ajout de vectux, vectuy, vectvy, vectvz double precision deterA, deltve ! ajout de deterA, deltve double precision epsrgp, blencp, climgp, epsilp, extrap, epsrsp double precision aa , aaa , aaaa , relaxp, thetap !=============================================================================== !=============================================================================== ! 0. GESTION MEMOIRE !=============================================================================== idebia = idbia0 idebra = idbra0 !=============================================================================== ! 1. PARAMETRAGE DU SOLVEUR ET INITIALISATION !=============================================================================== !--> Gradient Conjugue ireso1 = 0 !--> Parametrage de CODITS ! IVAR0= 0 LA VARIABLE N'EST ICI NI RIJ NI VITESSE ivar0 = 0 nitmap = 1000 ! IMRGRA = 0 nswrsp = 1 nswrgp = 100 imligp = -1 ircflp = 1 ischcp = 1 isstpp = 0 iescap = 0 imgr1 = 0 ncymap = 100 nitmgp = 10 iwarnp = iimlum blencp = zero epsilp = 1.d-8 epsrsp = 1.d-8 epsrgp = 1.d-5 climgp = 1.5d0 extrap = zero relaxp = 1.d0 !--> Il y a des dirichlets ndirc1 = 1 !--> Pas de convection pour le modele P1 iconv1 = 0 !--> Equation de diffusion idiff1 = 1 !--> Remise a zero des tableaux avant resolution do iel = 1,ncel drtp(iel) = zero theta4(iel) = zero thetaa(iel) = zero enddo do ifac = 1,nfac flurds(ifac) = zero enddo do ifac = 1,nfabor flurdb(ifac) = zero enddo !=============================================================================== ! 2. COEFFICIENT DE DIFFUSION AUX FACES !=============================================================================== do iel = 1,ncel ckmel(iel) = 1.d0 / ckmel(iel) enddo call viscfa & !========== ( idebia , idebra , & ndim , ncelet , ncel , nfac , nfabor , nfml , & nprfml , nnod , lndfac , lndfbr , ncelbr , & nideve , nrdeve , nituse , nrtuse , imvisf , & ifacel , ifabor , ifmfbr , ifmcel , iprfml , & ipnfac , nodfac , ipnfbr , nodfbr , & idevel , ituser , ia , & xyzcen , surfac , surfbo , cdgfac , cdgfbo , xyznod , & volume , ckmel , viscf , viscb , & rdevel , rtuser , ra ) !=============================================================================== ! 3. RESOLUTION !=============================================================================== ! Parametre pour schemas en temps et stationnaire thetap = 1.d0 idtva0 = 0 CNOM = ' ' WRITE(CNOM,'(A)') 'Rayon P1' inum = 1 nomvar(inum) = cnom call codits & !========== ( idebia , idebra , & ndim , ncelet , ncel , nfac , nfabor , nfml , nprfml , & nnod , lndfac , lndfbr , ncelbr , & nvar , nscal , nphas , & nideve , nrdeve , nituse , nrtuse , & idtva0 , ivar0 , iconv1 , idiff1 , ireso1 , ndirc1 , nitmap , & imrgra , nswrsp , nswrgp , imligp , ircflp , & ischcp , isstpp , iescap , & imgr1 , ncymap , nitmgp , inum , iwarnp , & blencp , epsilp , epsrsp , epsrgp , climgp , extrap , & relaxp , thetap , & ifacel , ifabor , ifmfbr , ifmcel , iprfml , & ipnfac , nodfac , ipnfbr , nodfbr , & idevel , ituser , ia , & xyzcen , surfac , surfbo , cdgfac , cdgfbo , xyznod , & volume , & thetaa , thetaa , cofrua , cofrub , cofrua , cofrub , & flurds , flurdb , & viscf , viscb , viscf , viscb , & rovsdt , smbrs , theta4 , & dam , xam , drtp , & w1 , w2 , w3 , w4 , w5 , & w6 , w7 , w8 , w9 , & rdevel , rtuser , ra ) !=============================================================================== ! 4. Vecteur densite de flux radiatif !=============================================================================== ! En periodique et parallele, echange avant calcul du gradient ! Parallele if (irangp.ge.0) then call parcom (theta4) !========== endif ! Periodique if (iperio.eq.1) then idimte = 0 itenso = 0 call percom & !========== ( idimte , itenso , & theta4 , theta4 , theta4 , & theta4 , theta4 , theta4 , & theta4 , theta4 , theta4) endif ! Calcul de la densite du flux radiatif QX, QY, QZ inc = 1 iccocg = 1 imligp = -1 iwarnp = iimlum epsrgp = 1.d-8 climgp = 1.5d0 extrap = 0.d0 nswrgp = 100 ivar0 = 0 iphydp = 0 call grdcel & !========== ( idebia , idebra , & ndim , ncelet , ncel , nfac , nfabor , nfml, & nprfml , & nnod , lndfac , lndfbr , ncelbr , nphas , & nideve , nrdeve , nituse , nrtuse , & ivar0 , imrgra , inc , iccocg , nswrgp , imligp, & iphydp , & iwarnp , nfecra , epsrgp , climgp , extrap , & ifacel , ifabor , ifmfbr , ifmcel , iprfml , & ipnfac , nodfac , ipnfbr , nodfbr , & idevel , ituser , ia , & xyzcen , surfac , surfbo , cdgfac , cdgfbo , xyznod , & volume , & w7 , w7 , w7 , & theta4 , cofrua , cofrub , & w1 , w2 , w3 , & w4 , w5 , w6 , & rdevel , rtuser , ra ) aa = - stephn * 4.d0 / 3.d0 do iel = 1,ncel aaa = aa * ckmel(iel) qx(iel) = w1(iel) * aaa qy(iel) = w2(iel) * aaa qz(iel) = w3(iel) * aaa enddo !=============================================================================== ! 5. Terme Source Radiatif d'absorption et densite de flux incident !=============================================================================== ! Calcul de la part d'absorption du terme Source Radiatif aa = 4.d0 * stephn do iel = 1,ncel sa(iel) = aa * theta4(iel) enddo ! Calcul du flux incident Qincid do ifac = 1, nfabor iel = ifabor(ifac) if (itypfb(ifac,iphas).eq.iparoi .or. & itypfb(ifac,iphas).eq.iparug ) then !--> Premiere version plus chere et legerement plus precise aaaa = tparoi(ifac)**4 aaa = 1.5d0 * ra(idistb-1+ifac) / ckmel(iel) & * ( 2.d0 /(2.d0-eps(ifac)) -1.d0 ) aa = ( aaa * aaaa + theta4(iel) ) / (1.d0 + aaa) qincid(ifac) = stephn * (2.d0 * aa - eps(ifac) * aaaa) & / (2.d0 - eps(ifac)) !--> Deuxieme version plus cheap mais moins precise ! QINCID(IFAC) = STEPHN * ! & (2.D0 * THETA4(IFABOR(IFAC)) - EPS(IFAC) * TPAROI(IFAC)**4) ! & / (2.D0 - EPS(IFAC)) else qincid(ifac) = stephn * theta4(iel) & + ( qx(iel) * surfbo(1,ifac) + & qy(iel) * surfbo(2,ifac) + & qz(iel) * surfbo(3,ifac) ) / & (0.5d0 * ra(isrfbn-1+ifac) ) endif enddo !=============================================================================== !6. Improved Differential Approximation (IDA) !=============================================================================== ! This module will improve the accuracy of P1-model near the wall and when the medium is optically thin. do iel = 1,ncel do fel = 1,nfabor !test pour voir si une paroi fait obstacle ou non -> détermination de visible : 1 veut dire que la paroi voit l'élément et 0 veut dire que la paroi ne voit pas l'élément. ! calcul du coefficient Delta qui représente la longueur du vecteur MI (M point du milieu et I point intercection entre MP et plan paroi). Si Delta est compris entre 0 et 1, cela veut dire qu'une paroi fait intercection à MP. Si Delta est supérieur à 1 ou inférieur à 0, cela veut dire que la paroi n'est pas entre M et P (elle est avant ou après). xpt1 = xyzcen(1,iel) ! x du point central de l'élément ypt1 = xyzcen(2,iel) ! y du point central de l'élément zpt1 = xyzcen(3,iel) ! z du point central de l'élément xpt2 = cdgfbo(1,fel) ! x du point central de la paroi ypt2 = cdgfbo(2,fel) ! y du point central de la paroi zpt2 = cdgfbo(3,fel) ! z du point central de la paroi visible = 1 ! Au départ, on considère que la paroi voit l'élément, cela sera modifié si une paroi fait obstacle numpar = 1 ! numéro de la première paroi testée en obstacle do while (numpar.le.nfabor .and. visible.eq.1) ! on s'arrête uniquement si une paroi est détectée en obstacle (visible = 0) ou si on a testé toutes les parois (numpar > nfabor) if (numpar.ne.fel) then ! on ne teste pas en obstacle la paroi considérée au départ xpt3 = cdgfbo(1,numpar) ! x du point central de la paroi d'intercection ypt3 = cdgfbo(2,numpar) ! y du point central de la paroi d'intercection zpt3 = cdgfbo(3,numpar) ! z du point central de la paroi d'intercection if (xpt3.ge.min(xpt1,xpt2) .and. xpt3.le.max(xpt1,xpt2) & ! test pour éviter de chercher les parois qui ne peuvent pas intercepter, si la position du milieu de la paroi n'est pas compris entre les positions de l'élement du milieu et de la paroi, ça ne peut pas être une paroi qui fait obstacle .and. ypt3.ge.min(ypt1,ypt2) .and. ypt3.le.max(ypt1,ypt2) & .and. zpt3.ge.min(zpt1,zpt2) .and. zpt3.le.max(zpt1,zpt2)) then vectux = surfbo(1,numpar) ! vecteur u de la surface (composante x) = vecteur x vectuy = surfbo(2,numpar)/2 ! vecteur u de la surface (composante y) = 1/2 vecteur y vectvy = surfbo(2,numpar)/2 ! vecteur v de la surface (composante y) = 1/2 vecteur y vectvz = surfbo(3,numpar) ! vecteur v de la surface (composante z) = vecteur z ! Le vecteur u n'a pas de composante en z et le vecteur v n'a pas de composante en x d'où l'absence de vectuz et vectvx deterA = (xpt2-xpt1)*vectuy*vectvz & ! déterminant de la matrice A -(ypt2-ypt1)*vectux*vectvz & +(zpt2-zpt1)*vectux*vectvy if (deterA.ne.0.0d0) then ! test au cas où le determinant est nul (paroi parallèle à la droite), si le déterminant n'est pas nul, on calcule le Delta Deltve = 1/deterA*(vectuy*vectvz*(xpt3-xpt1) & ! Calcul du Delta -vectux*vectvz*(ypt3-ypt1) & +vectux*vectvy*(zpt3-zpt1)) if (Deltve.ge.0.0d0 .and. Deltve.le.1.0d0) then ! Si le delta est entre 0 et 1, c'est que la paroi fait obstacle visible=0 ! visible = 0 veut dire obstacle endif endif endif endif numpar = numpar + 1 ! on incrémente le numéro de la paroi testée en obstacle enddo ! fin de boucle sur toutes les parois ou fin de boucle parce qu'une paroi fait obstacle ! Le reste n'est pas encore développer donc à ne pas prendre en considération pour le moment !if (visible.eq.1) then ! si visible vaut 1, cela veut dire que la paroi voit l'élément, on peut donc faire le calcul IDA. ! if (itypfb(fel,iphas).eq.isymet) then ! si la face est une symétrie !seeksymetry(xyzsym,nfacsy) ! ressort les coordonnées réelles du point symétrique et le numero de la face symetrique !!L = sqrt((xyzsym(1)-xyzcen(1,iel))**2+(xyzsym(2)-xyzcen(2,iel))**2+(xyzsym(3)-xyzcen(3,iel))**2) !longueur entre coordonnées réelles et élement ! else ! si la face n'est pas une symétrie (donc mur) voir si mettre un elseif(itypfb(fel,iphas).eq.iparoi .or. itypfb(fel,iphas).eq.iparug) !L = sqrt((cdgfbo(1,fel)-xyzcen(1,iel))**2+(cdgfbo(2,fel)-xyzcen(2,iel))**2+(cdgfbo(3,fel)-xyzcen(3,iel))**2) ! longueur entre paroi et élément !nfacsy = fel ! numero de la face symetrique = numero de la face ! endif !calcul cofext suivant le coef d extinction de chaque element traversé !tau = cofext*L ! coefficient d'extinction !s0 = (1.0d0-tau*exp(-tau))/(coefext*(1.0d0-exp(-tau))) ! longueur s0 !recherche élément à longueur s0(g,qx,qy,qz,numele) ! lorsque l'on a trouvé l'élément à une distance s0 dans la même direction, on récupère la valeur de la luminance et le flux radiatif en x,y,z dans cet élément !setoil = (1.0d0-propdf(numele)/cofext(numele))*eps(nfacsy)*stephn*tparoi(nfacsy)**4+propdf(numele)/(4*pi*cofext(numele))*(g+propfp(numele)*sqrt(qx**2+qy**2+qz**2)) !radios = eps(nfacsy)*stephn*tparoi(nfacsy)**4 ! endif ! fin du if sur visible, si visible = 0 , on ne fait pas le calcul IDA pour ce couple frontière-élément, on passe à une autre frontière. enddo ! fin de boucle sur toutes les frontières pour cet élément, on passe à l'élément suivant. enddo !=============================================================================== !-------- ! FORMATS !-------- !---- ! FIN !---- return end subroutine