daasolvp.F

Go to the documentation of this file.

Copyright>        OpenRadioss
Copyright>        Copyright (C) 1986-2025 Altair Engineering Inc.
Copyright>
Copyright>        This program is free software: you can redistribute it and/or modify
Copyright>        it under the terms of the GNU Affero General Public License as published by
Copyright>        the Free Software Foundation, either version 3 of the License, or
Copyright>        (at your option) any later version.
Copyright>
Copyright>        This program is distributed in the hope that it will be useful,
Copyright>        but WITHOUT ANY WARRANTY; without even the implied warranty of
Copyright>        MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
Copyright>        GNU Affero General Public License for more details.
Copyright>
Copyright>        You should have received a copy of the GNU Affero General Public License
Copyright>        along with this program.  If not, see <https://www.gnu.org/licenses/>.
Copyright>
Copyright>
Copyright>        Commercial Alternative: Altair Radioss Software
Copyright>
Copyright>        As an alternative to this open-source version, Altair also offers Altair Radioss
Copyright>        software under a commercial license.  Contact Altair to discuss further if the
Copyright>        commercial version may interest you: https://www.altair.com/radioss/.
!||====================================================================
!||    daasolvp               ../engine/source/fluid/daasolvp.F
!||--- called by ------------------------------------------------------
!||    flow0                  ../engine/source/fluid/flow0.F
!||--- calls      -----------------------------------------------------
!||    finter                 ../engine/source/tools/curve/finter.F
!||    finter_smooth          ../engine/source/tools/curve/finter_smooth.f
!||    spmd_fl_sum            ../engine/source/mpi/generic/spmd_fl_sum.F
!||--- uses       -----------------------------------------------------
!||    python_funct_mod       ../common_source/modules/python_mod.f90
!||    sensor_mod             ../common_source/modules/sensor_mod.F90
!||====================================================================
      SUBROUTINE daasolvp(NDIM,   NNO,    NEL,   NEL_LOC,
     .                    IFLOW,  IBUF,   ELEM,  IBUFL,   SHELL_GA,  CNP,
     .                    RFLOW,  NORMAL, TA,    AREAF,   COSG,      DIST,  
     .                    MFLE,   ACCF,   PM,    PTI, X,  V, A, NPC, TF, 
     .                    NBGAUGE,LGAUGE, GAUGE, NSENSOR, SENSOR_TAB, IGRV, AGRV,
     .                    CBEM   ,IPIV_L, MFLE_L,IBUFELR, IBUFELC,NFUNCT,
     .                    PYTHON ,WFEXT)
C-----------------------------------------------
C   M o d u l e s
C-----------------------------------------------  
      USE sensor_mod
      USE python_funct_mod
C-----------------------------------------------
C   I m p l i c i t   T y p e s
C-----------------------------------------------
#include      "implicit_f.inc"
C-----------------------------------------------
C   C o m m o n   B l o c k s
C-----------------------------------------------
#include      "param_c.inc"
!#include      "com06_c.inc"
#include      "com08_c.inc"
C-----------------------------------------------
C   D u m m y   A r g u m e n t s
C-----------------------------------------------
      INTEGER NDIM, NNO, NEL, NBGAUGE, NSENSOR, NEL_LOC,NFUNCT
      INTEGER IFLOW(*), IBUF(*), ELEM(NDIM,*), SHELL_GA(*), NPC(*), LGAUGE(3,*)
      INTEGER IBUFL(*), CNP(*), IGRV(NIGRV,*)
      INTEGER IPIV_L(*), IBUFELR(NEL), IBUFELC(NEL)
      my_real X(3,*), V(3,*), A(3,*), TF(*), RFLOW(*)
      my_real NORMAL(3,*), TA(*), AREAF(*), COSG(*), DIST(*), PM(*), ACCF(*), PTI(*)
      my_real mfle(nel,*), gauge(llgauge,*), agrv(lfacgrv,*)
      my_real cbem(nel,*), mfle_l(nel_loc,*)
      TYPE (SENSOR_STR_) ,DIMENSION(NSENSOR) :: SENSOR_TAB
      TYPE (PYTHON_), INTENT(INOUT) :: PYTHON
      DOUBLE PRECISION,INTENT(INOUT) :: WFEXT
C-----------------------------------------------
C   L o c a l   V a r i a b l e s
C-----------------------------------------------
      INTEGER  I, J, K, L, I1, I2, I3, I4, I5, N1, N2, N3, N4, II, JJ, KK
      INTEGER  IFPRES, JFORM, KFORM, IPRES, IWAVE, INTEGR, FREESURF, AFTERFLOW
      INTEGER  GRAV_ID, IFUNC, ISENS, LENBUF, NNO_L, NEL_L, NELMAX, INFO,ISMOOTH
      INTEGER  NEBUF(NEL)
      INTEGER  NBLOC, NPROW, NPCOL, ICTXT, DESCA(9), DESCB(9), MYROW, MYCOL, MXLLDM, NUMROC
      EXTERNAL NUMROC
      my_real  x1,  y1,  z1,  x2,  y2,  z2,  x3,  y3,  z3,  x4,  y4,  z4,
     .         x0,  y0,  z0,  x12, y12, z12, x13, y13, z13, x24, y24, z24,    
     .         nrx, nry, nrz, area2, wf(2), wi(4,2),
     .         pm1, dsc, vx0, vy0, vz0
      my_real  rho, ssp, rhoc, tta, dydx, ff, sfpres, dppdt, pmax, theta, norm, fac1, ts, gravity
      my_real  dirwix, dirwiy, dirwiz, dirwrx, dirwry, dirwrz
      my_real  off(nel), pp(nel), vi(nel), ps(nel), fel(nel), pr(nel), vr(nel)
      my_real  dpidt(nel), dprdt(nel), dpmdt(nel), hh(nel)
      my_real  xa,  ya,  za,  fcx,  fcy, pmin, ptot
      my_real  apmax, atheta, ratio
      my_real  finter,finter_smooth
      EXTERNAL finter,finter_smooth
      my_real  vect(nel), vect1(nel), vel(nel)
      my_real  xl(3,nno), vl(3,nno)
      my_real, ALLOCATABLE :: sbuf(:), rbuf(:)
 
#if defined(MPI) && defined(MYREAL8) && !defined(WITHOUT_LINALG)
C Coordonnees et vitesses locales SPMD
      nno_l=iflow(16)
      lenbuf=6*nno
      ALLOCATE(sbuf(lenbuf), rbuf(lenbuf))
      sbuf(1:lenbuf)=zero
      rbuf(1:lenbuf)=zero
      DO i=1,nno_l
         ii=ibufl(i)
         jj=ibuf(ii)
         kk=6*(ii-1)
         sbuf(kk+1)=x(1,jj)/cnp(ii)
         sbuf(kk+2)=x(2,jj)/cnp(ii)
         sbuf(kk+3)=x(3,jj)/cnp(ii)
         sbuf(kk+4)=v(1,jj)/cnp(ii)
         sbuf(kk+5)=v(2,jj)/cnp(ii)
         sbuf(kk+6)=v(3,jj)/cnp(ii)
      ENDDO
 
      CALL spmd_fl_sum(sbuf, lenbuf, rbuf)
 
      DO i=1,nno
         k=6*(i-1)
         xl(1,i)=rbuf(k+1)
         xl(2,i)=rbuf(k+2)
         xl(3,i)=rbuf(k+3)
         vl(1,i)=rbuf(k+4)
         vl(2,i)=rbuf(k+5)
         vl(3,i)=rbuf(k+6)
      ENDDO
      DEALLOCATE(sbuf, rbuf)
 
      jform  = iflow(4)
      nbloc  = iflow(12)
      nprow  = iflow(18)
      npcol  = iflow(19)
      ipres  = iflow(21)
      iwave  = iflow(22)
      kform  = iflow(23)
      integr = iflow(24)
      freesurf  = iflow(25)
      afterflow = iflow(26)
      grav_id   = iflow(27)
      nelmax    = iflow(28)
      rhoc   = rflow(1)
      ssp    = rflow(2)
      rho    = rflow(5)
      fac1   = rflow(8)
      dsc    = rflow(12)
      xa     = rflow(16)
      ya     = rflow(17)
      za     = rflow(18)
      pmin   = rflow(22)
      apmax  = rflow(23)
      atheta = rflow(24)
 
      IF(jform == 1) THEN
        DO i = 1,nel
           i1 = elem(1,i)
           i2 = elem(2,i)
           i3 = elem(3,i)
           x1 = xl(1,i1)
           x2 = xl(1,i2)
           x3 = xl(1,i3)
           y1 = xl(2,i1)
           y2 = xl(2,i2)
           y3 = xl(2,i3)
           z1 = xl(3,i1)
           z2 = xl(3,i2)
           z3 = xl(3,i3)
           x12= x2-x1
           y12= y2-y1
           z12= z2-z1
           x13= x3-x1
           y13= y3-y1
           z13= z3-z1
           nrx= y12*z13-z12*y13
           nry= z12*x13-x12*z13
           nrz= x12*y13-y12*x13
           area2 = sqrt(nrx**2+nry**2+nrz**2)
           normal(1,i) = nrx/area2
           normal(2,i) = nry/area2
           normal(3,i) = nrz/area2
           areaf(i) = half*area2
           IF(iwave == 1)THEN
C Spherical wave
             x0=third*x1+third*x2+third*x3
             y0=third*y1+third*y2+third*y3
             z0=third*z1+third*z2+third*z3
             dirwix=x0-rflow(9)
             dirwiy=y0-rflow(10)
             dirwiz=z0-rflow(11)
             norm=sqrt(dirwix**2+dirwiy**2+dirwiz**2)
             dirwix=dirwix/norm
             dirwiy=dirwiy/norm
             dirwiz=dirwiz/norm
             IF(freesurf == 2) THEN
               dirwrx=x0-rflow(13)
               dirwry=y0-rflow(14)
               dirwrz=z0-rflow(15)
               norm=sqrt(dirwrx**2+dirwry**2+dirwrz**2)
               dirwrx=dirwrx/norm
               dirwry=dirwry/norm
               dirwrz=dirwrz/norm
                     cosg(i+nel)=normal(1,i)*dirwrx+normal(2,i)*dirwry+normal(3,i)*dirwrz
             ENDIF
C hydrostatic pressure
             hh(i)=zero
             IF(grav_id > 0) THEN
               hh(i)=max(zero,(xa-x0)*rflow(19)+(ya-y0)*rflow(20)+(za-z0)*rflow(21))
             ENDIF
           ELSEIF(iwave == 2) THEN
C Onde plane
             dirwix=rflow(9)
             dirwiy=rflow(10)
             dirwiz=rflow(11)
C hydrostatic pressure
             hh(i)=zero
           ENDIF
                 cosg(i)=normal(1,i)*dirwix+normal(2,i)*dirwiy+normal(3,i)*dirwiz
C vitesse N-1/2
           vx0 = third*vl(1,i1)+third*vl(1,i2)+third*vl(1,i3)
           vy0 = third*vl(2,i1)+third*vl(2,i2)+third*vl(2,i3)
           vz0 = third*vl(3,i1)+third*vl(3,i2)+third*vl(3,i3)
           vel(i) = vx0*normal(1,i)+vy0*normal(2,i)+vz0*normal(3,i)
        ENDDO
      ELSEIF(jform == 2) THEN
        wi(1,1)=fourth
        wi(2,1)=fourth
        wi(3,1)=fourth
        wi(4,1)=fourth
        wi(1,2)=third
        wi(2,2)=third
        wi(3,2)=one_over_6
        wi(4,2)=one_over_6
        DO i = 1,nel
           i1 = elem(1,i)
           i2 = elem(2,i)
           i3 = elem(3,i)
           i4 = elem(4,i)
           i5 = elem(5,i)
           x1 = xl(1,i1)
           x2 = xl(1,i2)
           x3 = xl(1,i3)
           x4 = xl(1,i4)
           y1 = xl(2,i1)
           y2 = xl(2,i2)
           y3 = xl(2,i3)
           y4 = xl(2,i4)
           z1 = xl(3,i1)
           z2 = xl(3,i2)
           z3 = xl(3,i3)
           z4 = xl(3,i4)
           x13= x3-x1
           y13= y3-y1
           z13= z3-z1
           x24= x4-x2
           y24= y4-y2
           z24= z4-z2
           nrx=y13*z24-z13*y24
           nry=z13*x24-x13*z24
           nrz=x13*y24-y13*x24
           area2 = sqrt(nrx**2+nry**2+nrz**2)
           normal(1,i) = nrx/area2
           normal(2,i) = nry/area2
           normal(3,i) = nrz/area2
           areaf(i) = half*area2
           IF(iwave == 1)THEN
C Spherical wave
             x0=wi(1,i5)*x1+wi(2,i5)*x2+wi(3,i5)*x3+wi(4,i5)*x4
             y0=wi(1,i5)*y1+wi(2,i5)*y2+wi(3,i5)*y3+wi(4,i5)*y4
             z0=wi(1,i5)*z1+wi(2,i5)*z2+wi(3,i5)*z3+wi(4,i5)*z4
             dirwix=x0-rflow(9)
             dirwiy=y0-rflow(10)
             dirwiz=z0-rflow(11)
             norm=sqrt(dirwix**2+dirwiy**2+dirwiz**2)
             dirwix=dirwix/norm
             dirwiy=dirwiy/norm
             dirwiz=dirwiz/norm
             IF(freesurf == 2) THEN
               dirwrx=x0-rflow(13)
               dirwry=y0-rflow(14)
               dirwrz=z0-rflow(15)
               norm=sqrt(dirwrx**2+dirwry**2+dirwrz**2)
               dirwrx=dirwrx/norm
               dirwry=dirwry/norm
               dirwrz=dirwrz/norm
                     cosg(i+nel)=normal(1,i)*dirwrx+normal(2,i)*dirwry+normal(3,i)*dirwrz
             ENDIF
C hydrostatic pressure
             hh(i)=zero
             IF(grav_id > 0) THEN
               hh(i)=max(zero,(xa-x0)*rflow(19)+(ya-y0)*rflow(20)+(za-z0)*rflow(21))
             ENDIF
           ELSEIF(iwave == 2) THEN
C Onde plane
             dirwix=rflow(9)
             dirwiy=rflow(10)
             dirwiz=rflow(11)
C hydrostatic pressure
             hh(i)=zero
           ENDIF
           cosg(i)=normal(1,i)*dirwix+normal(2,i)*dirwiy+normal(3,i)*dirwiz
C vitesse N-1/2
           vx0 = wi(1,i5)*vl(1,i1)+wi(2,i5)*vl(1,i2)+wi(3,i5)*vl(1,i3)+wi(4,i5)*vl(1,i4)
           vy0 = wi(1,i5)*vl(2,i1)+wi(2,i5)*vl(2,i2)+wi(3,i5)*vl(2,i3)+wi(4,i5)*vl(2,i4)
           vz0 = wi(1,i5)*vl(3,i1)+wi(2,i5)*vl(3,i2)+wi(3,i5)*vl(3,i3)+wi(4,i5)*vl(3,i4)
           vel(i) = vx0*normal(1,i)+vy0*normal(2,i)+vz0*normal(3,i)
        ENDDO
 
      ENDIF
C-----------------------------------------------
C    1. Pression incidente
C-----------------------------------------------
      DO i=1,nel
         off(i) = zero
         tta = tt-ta(i)
         IF(tta > zero) off(i) = one
      ENDDO
      IF(ipres == 1) THEN
        IF(iwave == 1)THEN
          DO i=1,nel
             tta   = tt-ta(i)
             ratio = dsc/dist(i)
             pmax  = rflow(6)*ratio**apmax
             theta = rflow(7)*ratio**atheta
             pp(i) = pmax*exp(-tta/theta)*off(i)
             dpidt(i) = -pp(i)*cosg(i)/theta
          ENDDO
        ELSEIF(iwave == 2) THEN
          pmax  = rflow(6)
          theta = rflow(7)
          DO i=1,nel
             tta = tt-ta(i)
             pp(i) = pmax*exp(-tta/theta)*off(i)
             dpidt(i) = -pp(i)*cosg(i)/theta
          ENDDO
        ENDIF
      ELSEIF(ipres == 2) THEN
        ifpres = iflow(7)
        sfpres = rflow(3)
        IF(iwave == 1)THEN
          DO i=1,nel
             tta   = tt-ta(i)
             ratio = dsc/dist(i)
             IF(ifpres > 0) THEN
               ismooth = npc(2*nfunct+ifpres+1)
!!                PP(I) = SFPRES*FINTER(IFPRES,TTA,NPC,TF,DPPDT)*RATIO*OFF(I)
                IF (ismooth == 0) THEN
                  pp(i) = sfpres*finter(ifpres,tta,npc,tf,dppdt)*ratio*off(i)
                ELSE IF(ismooth > 0) THEN
                  pp(i) = sfpres*finter_smooth(ifpres,tta,npc,tf,dppdt)*ratio*off(i)
                ELSE 
                  ismooth = -ismooth ! function id
                  CALL python_call_funct1d(python, ismooth,tta,pp(i))
                  pp(i) = pp(i)*sfpres*ratio*off(i)
                  CALL python_deriv_funct1d(python, ismooth,tta,dppdt)
                ENDIF ! IF (ISMOOTH == 0)
                dpidt(i) = dppdt*cosg(i)*ratio*off(i)
             ELSE
                pp(i) = sfpres*ratio*off(i)
                dpidt(i) = zero
             ENDIF
          ENDDO
        ELSEIF(iwave == 2) THEN
          DO i=1,nel
             tta = tt-ta(i)
             IF(ifpres > 0) THEN
               ismooth = npc(2*nfunct+ifpres+1)
!!                PP(I) = SFPRES*FINTER(IFPRES,TTA,NPC,TF,DPPDT)*OFF(I)
                IF (ismooth == 0) THEN
                  pp(i) = sfpres*finter(ifpres,tta,npc,tf,dppdt)*off(i)
                ELSE IF (ismooth > 0) THEN
                  pp(i) = sfpres*finter_smooth(ifpres,tta,npc,tf,dppdt)*off(i)
                ELSE 
                  ismooth = -ismooth ! function id
                  CALL python_call_funct1d(python, ismooth,tta,pp(i))
                  pp(i) = pp(i)*sfpres*off(i)
                  CALL python_deriv_funct1d(python, ismooth,tta,dppdt)
                ENDIF ! IF (ISMOOTH == 0)
                dpidt(i) = dppdt*cosg(i)*off(i)
             ELSE
                pp(i) = sfpres*off(i)
                dpidt(i) = zero
             ENDIF
          ENDDO
        ENDIF
      ENDIF    
C-----------------------------------------------
C     Vitesse onde incidente
C-----------------------------------------------
      DO i=1,nel
         vi(i)=pp(i)*cosg(i)/rhoc
      ENDDO
      IF(afterflow == 2) THEN
         IF(iwave == 1)THEN
C add afterflow velocity, no afterflow velocity for plane wave
           IF(ipres == 1) THEN
              DO i=1,nel
                 ratio = dsc/dist(i)
                 pmax  = rflow(6)*ratio**apmax
                 theta = rflow(7)*ratio**atheta
                 vi(i) = vi(i) + (pmax-pp(i))*cosg(i)*theta/(rho*dist(i))
              ENDDO
           ELSEIF(ipres == 2) THEN
              DO i=1,nel
                 pti(i) = pti(i) + pp(i)*dt1
                   vi(i) = vi(i) + pti(i)*cosg(i)/(rho*dist(i))
              ENDDO
           ENDIF
         ENDIF
      ENDIF
C-----------------------------------------------
C   2. Pression Onde Reflechie
C-----------------------------------------------
      IF(freesurf == 2) THEN
        DO i=1,nel
           off(i) = zero
           tta = tt-ta(i+nel)
           IF(tta > zero) off(i) = one
        ENDDO
        IF(ipres == 1) THEN
          DO i=1,nel
             j = i+nel
             tta = tt-ta(j)
             ratio = dsc/dist(j)
             pmax  = rflow(6)*ratio**apmax
             theta = rflow(7)*ratio**atheta
             pr(i) = -pmax*exp(-tta/theta)*off(i)
             dprdt(i) = -pr(i)*cosg(j)/theta
          ENDDO
        ELSEIF(ipres == 2) THEN
          DO i=1,nel
             j = i+nel
             tta = tt-ta(j)
             ratio = dsc/dist(j)
             IF(ifpres > 0) THEN
               ismooth = npc(2*nfunct+ifpres+1)
!!                PR(I) = -SFPRES*FINTER(IFPRES,TTA,NPC,TF,DPPDT)*RATIO*OFF(I)
                IF (ismooth == 0) THEN
                  pr(i) = -sfpres*finter(ifpres,tta,npc,tf,dppdt)*ratio*off(i)
                ELSE IF (ismooth > 0) THEN
                  pr(i) = -sfpres*finter_smooth(ifpres,tta,npc,tf,dppdt)*ratio*off(i)
                ELSE 
                  ismooth = -ismooth ! function id
                  CALL python_call_funct1d(python, ismooth,tta,pr(i))
                  pr(i) = -pr(i)*sfpres*ratio*off(i)
                  CALL python_deriv_funct1d(python, ismooth,tta,dppdt)
                ENDIF ! IF (ISMOOTH == 0)
                dprdt(i) = dppdt*cosg(j)*off(i)
             ELSE
                pr(i) = -sfpres*ratio*off(i)
                dprdt(i) = zero
             ENDIF
          ENDDO
        ENDIF    
C-----------------------------------------------
C    Vitesse onde reflechie
C-----------------------------------------------
        DO i=1,nel
           vr(i)=pr(i)*cosg(i+nel)/rhoc
        ENDDO
      ELSE  ! pas de surface libre
        DO i=1,nel
          pr(i) = zero
          vr(i) = zero
          dprdt(i) = zero
        ENDDO
      ENDIF
C-----------------------------------------------
C    3. Scattered pressure Schema Euler ordre 1
C-----------------------------------------------
      IF(kform == 1) THEN
C DAA-1
        IF(nel < nelmax) THEN
          DO i=1,nel
             dpmdt(i) = rhoc*accf(i)
             vect(i)  = rhoc*areaf(i)*(pm(i)-rhoc*(vi(i)+vr(i)))
          ENDDO
          vect(1:nel) = matmul(mfle(1:nel, 1:nel), vect(1:nel))
          DO i=1,nel
             dpmdt(i) = dpmdt(i) - vect(i)
          ENDDO
          IF(integr == 1) THEN
            DO i=1,nel
               pm(i) = pm(i) + dpmdt(i)*dt1
               ps(i) = pm(i) - rhoc*(vi(i)+vr(i))
            ENDDO
 
          ELSEIF(integr == 2) THEN
C Prediction
            DO i=1,nel
               ps(i) = pm(i) + dpmdt(i)*dt1*half - rhoc*(vi(i)+vr(i))
            ENDDO
C Correction
            DO i=1,nel
               dpmdt(i) = rhoc*accf(i)
               vect(i)  = rhoc*areaf(i)*ps(i)
            ENDDO
            vect(1:nel) = matmul(mfle(1:nel, 1:nel), vect(1:nel))
            DO i=1,nel
               dpmdt(i) = dpmdt(i) - vect(i)
            ENDDO
            DO i=1,nel
               pm(i) = pm(i) + dpmdt(i)*dt1
               ps(i) = pm(i) - rhoc*(vi(i)+vr(i)) 
            ENDDO
          ENDIF
C
        ELSE  ! resolution systeme  MFLE*X=Y
C  Initialisation de la process grid
          CALL sl_init(ictxt, nprow, npcol)
          CALL blacs_gridinfo(ictxt, nprow, npcol, myrow, mycol)
C          MXLLDM = NUMROC(NEL, NBLOC, MYROW, 1, NPROW)
          mxlldm = nel_loc
          CALL descinit(desca, nel, nel, nbloc, nbloc, 0, 0, ictxt, mxlldm, info)
          CALL descinit(descb, nel,   1, nbloc, nbloc, 0, 0, ictxt, mxlldm, info)
 
          IF(dt1 == zero) THEN
C  Compute local sub-matrix 
             k=0
             DO i=1,nel
                IF(ibufelr(i) == 0) cycle
                k=k+1
                l=0
                DO j=1,nel
                  IF(ibufelc(j) == 0) cycle
                  l=l+1
                  mfle_l(k,l) = mfle(i,j)
                ENDDO
             ENDDO
C  L U decomposition
             CALL pdgetrf(nel, nel, mfle_l, 1, 1, desca, ipiv_l, info)
          ENDIF
 
          DO i=1,nel
           vect(i) = two*ssp*(rhoc*(vi(i)+vr(i))-pm(i))
          ENDDO
          CALL dgemv('T',nel, nel, one, cbem, nel, vect, 1, zero, vect1, 1)
 
          vect(1:nel)=zero
          IF(ibufelc(1) > 0) THEN
            k=0
            DO i=1,nel
               IF(ibufelr(i) == 0) cycle
               k=k+1
               vect(k) = vect1(i)
            ENDDO
          ENDIF
            
C         Resolution
          CALL pdgetrs('N', nel, 1, mfle_l, 1, 1, desca, ipiv_l, vect, 1, 1, descb, info)
 
C Echange 
          ALLOCATE(sbuf(nel), rbuf(nel))
          sbuf(1:nel)=zero
          rbuf(1:nel)=zero
          IF(ibufelc(1) > 0) THEN
            k=0
            DO i=1,nel
               IF(ibufelr(i) == 0) cycle
               k=k+1
               sbuf(i) = vect(k)
            ENDDO
          ENDIF
          CALL spmd_fl_sum(sbuf, nel, rbuf)
          DO i=1,nel
             vect(i)=rbuf(i)
          ENDDO
 
          vect1(1:nel) = matmul(cbem(1:nel,1:nel), vect(1:nel))
          DO i=1,nel
             IF(pp(i) /= zero) THEN
               dpmdt(i) = rhoc*accf(i) + vect(i)
             ELSE
               dpmdt(i)=zero
             ENDIF
          ENDDO
 
          IF(integr == 1) THEN
             DEALLOCATE(sbuf, rbuf)
             DO i=1,nel
                pm(i) = pm(i) + dpmdt(i)*dt1
                ps(i) = pm(i) - rhoc*(vi(i)+vr(i))
             ENDDO
 
          ELSEIF(integr == 2) THEN
C Prediction
             DO i=1,nel
                ps(i) = pm(i) + dpmdt(i)*dt1*half - rhoc*(vi(i)+vr(i))
             ENDDO
C Correction
             DO i=1,nel
              vect(i) = -two*ssp*ps(i)
             ENDDO
             CALL dgemv('T',nel, nel, one, cbem, nel, vect, 1, zero, vect1, 1)
             vect(1:nel)=zero
             IF(ibufelc(1) > 0) THEN
               k=0
               DO i=1,nel
                  IF(ibufelr(i) == 0) cycle
                  k=k+1
                  vect(k) = vect1(i)
               ENDDO
             ENDIF
 
C Resolution
             CALL pdgetrs('N', nel, 1, mfle_l, 1, 1, desca, ipiv_l, vect, 1, 1, descb, info)
C Echange 
             sbuf(1:nel)=zero
             rbuf(1:nel)=zero
             IF(ibufelc(1) > 0) THEN
               k=0
               DO i=1,nel
                  IF(ibufelr(i) == 0) cycle
                  k=k+1
                  sbuf(i) = vect(k)
               ENDDO
             ENDIF
             CALL spmd_fl_sum(sbuf, nel, rbuf)
             DO i=1,nel
                vect(i)=rbuf(i)
             ENDDO
             DEALLOCATE(sbuf, rbuf)
 
             vect1(1:nel) = matmul(cbem(1:nel,1:nel), vect(1:nel))
             DO i=1,nel
                dpmdt(i) = rhoc*accf(i) + vect1(i)
             ENDDO
             DO i=1,nel
                pm(i) = pm(i) + dpmdt(i)*dt1
                ps(i) = pm(i) - rhoc*(vi(i)+vr(i))
             ENDDO
          ENDIF
 
          CALL blacs_gridexit(ictxt)
        ENDIF
C
      ELSEIF(kform == 2) THEN
C High Frequency Approximation
        DO i=1,nel
           dpmdt(i) = rhoc*accf(i)
        ENDDO
        DO i=1,nel
           pm(i) = pm(i) + dpmdt(i)*dt1
           ps(i) = pm(i) - rhoc*vi(i) 
        ENDDO
C
      ELSEIF(kform == 3) THEN
C Virtual Mass Approximation
        DO i=1,nel
           ps(i) = zero
           DO j=1,nel
              ps(i) = ps(i) + mfle(i,j)*(accf(j)-(dpidt(j)+dprdt(j))/rhoc)
           ENDDO
           ps(i) = ps(i)/areaf(i)
        ENDDO
      ENDIF
C-----------------------------------------------
C    4. gravity
C-----------------------------------------------
      gravity = zero
      IF(grav_id > 0) THEN
        fcy = agrv(1,grav_id)
        fcx = agrv(2,grav_id)
        ifunc = igrv(3,grav_id)
        isens = 0
        DO k=1,nsensor
          IF(igrv(6,grav_id) == sensor_tab(k)%SENS_ID) isens=k ! do it in starter !!!
        ENDDO
        IF(isens==0)THEN
          ts = tt
        ELSE
          ts = tt-sensor_tab(isens)%TSTART
        ENDIF
        ismooth = 0
        IF (ifunc > 0) THEN
          ismooth = npc(2*nfunct+ifunc+1)
!!        GRAVITY = FCY*FINTER(IFUNC,TS*FCX,NPC,TF,DYDX)
          IF (ismooth == 0) THEN
            gravity = fcy*finter(ifunc,ts*fcx,npc,tf,dydx)
          ELSE IF(ismooth > 0) THEN
            gravity = fcy*finter_smooth(ifunc,ts*fcx,npc,tf,dydx)
          ELSE 
            ismooth = -ismooth ! function id
            CALL python_call_funct1d(python, ismooth,tta,gravity)
            gravity = gravity*fcy
          ENDIF ! IF (ISMOOTH == 0)
        ELSE
          gravity   = fcy
        ENDIF
      ENDIF
C-----------------------------------------------
C    5. compute total pressure forces
C-----------------------------------------------
      DO i=1,nel
         ptot = pp(i)+ps(i)+pr(i)+rho*hh(i)*abs(gravity)
         fel(i) = areaf(i)*max(ptot,pmin)
      ENDDO
C    
      IF(jform == 1) THEN
        nel_l=0
        DO i=1,nel
           i1 = elem(1,i)
           i2 = elem(2,i)
           i3 = elem(3,i)
           n1 = ibuf(i1)
           n2 = ibuf(i2)
           n3 = ibuf(i3)
           IF(n1/=0.OR.n2/=0.OR.n3/=0) THEN
              nel_l = nel_l+1
              nebuf(nel_l) = i
           ENDIF
        ENDDO
        DO k=1,nel_l
           i  = nebuf(k)
           i1 = elem(1,i)
           i2 = elem(2,i)
           i3 = elem(3,i)
           n1 = ibuf(i1)
           n2 = ibuf(i2)
           n3 = ibuf(i3)
           nrx= normal(1,i)
           nry= normal(2,i)
           nrz= normal(3,i)
           ff = fac1*third*fel(i)
           IF(n1/=0) THEN
              a(1,n1) = a(1,n1) + ff*nrx
              a(2,n1) = a(2,n1) + ff*nry
              a(3,n1) = a(3,n1) + ff*nrz
              wfext   = wfext+ff*dt1*(nrx*v(1,n1)+nry*v(2,n1)+nrz*v(3,n1))/cnp(i1)
           ENDIF
           IF(n2/=0) THEN
              a(1,n2) = a(1,n2) + ff*nrx
              a(2,n2) = a(2,n2) + ff*nry
              a(3,n2) = a(3,n2) + ff*nrz
              wfext   = wfext+ff*dt1*(nrx*v(1,n2)+nry*v(2,n2)+nrz*v(3,n2))/cnp(i2)
           ENDIF
           IF(n3/=0) THEN
              a(1,n3) = a(1,n3) + ff*nrx
              a(2,n3) = a(2,n3) + ff*nry
              a(3,n3) = a(3,n3) + ff*nrz
              wfext   = wfext+ff*dt1*(nrx*v(1,n3)+nry*v(2,n3)+nrz*v(3,n3))/cnp(i3)
           ENDIF
        ENDDO
      ELSEIF(jform == 2) THEN
        nel_l=0
        DO i=1,nel
           i1 = elem(1,i)
           i2 = elem(2,i)
           i3 = elem(3,i)
           i4 = elem(4,i)
           n1 = ibuf(i1)
           n2 = ibuf(i2)
           n3 = ibuf(i3)
           n4 = ibuf(i4)
           IF(n1/=0.OR.n2/=0.OR.n3/=0.OR.n4/=0) THEN
              nel_l = nel_l+1
              nebuf(nel_l) = i
           ENDIF
        ENDDO
        wf(1)=fac1*fourth
        wf(2)=fac1*third
        DO k=1,nel_l
           i  = nebuf(k)
           i1 = elem(1,i)
           i2 = elem(2,i)
           i3 = elem(3,i)
           i4 = elem(4,i)
           i5 = elem(5,i)
           n1 = ibuf(i1)
           n2 = ibuf(i2)
           n3 = ibuf(i3)
           n4 = ibuf(i4)
           nrx= normal(1,i)
           nry= normal(2,i)
           nrz= normal(3,i)
           ff = wf(i5)*fel(i)
           IF(n1/=0) THEN
              a(1,n1) = a(1,n1) + ff*nrx
              a(2,n1) = a(2,n1) + ff*nry
              a(3,n1) = a(3,n1) + ff*nrz
              wfext   = wfext+ff*dt1*(nrx*v(1,n1)+nry*v(2,n1)+nrz*v(3,n1))/cnp(i1)
           ENDIF
           IF(n2/=0) THEN
              a(1,n2) = a(1,n2) + ff*nrx
              a(2,n2) = a(2,n2) + ff*nry
              a(3,n2) = a(3,n2) + ff*nrz
              wfext   = wfext+ff*dt1*(nrx*v(1,n2)+nry*v(2,n2)+nrz*v(3,n2))/cnp(i2)
           ENDIF
           IF(n3/=0) THEN
              a(1,n3) = a(1,n3) + ff*nrx
              a(2,n3) = a(2,n3) + ff*nry
              a(3,n3) = a(3,n3) + ff*nrz
              wfext   = wfext+ff*dt1*(nrx*v(1,n3)+nry*v(2,n3)+nrz*v(3,n3))/cnp(i3)
           ENDIF
           IF(n4/=0.AND.i5==1) THEN
              a(1,n4) = a(1,n4) + ff*nrx
              a(2,n4) = a(2,n4) + ff*nry
              a(3,n4) = a(3,n4) + ff*nrz
              wfext   = wfext+ff*dt1*(nrx*v(1,n4)+nry*v(2,n4)+nrz*v(3,n4))/cnp(i4)
           ENDIF
        ENDDO
      ENDIF
C-----------------------------------------------------------
C     5. Gauge
C-----------------------------------------------------------
      IF(jform == 2) THEN
        DO i=1,nbgauge
           gauge(30,i)=zero
           IF(lgauge(1,i) == 0 .AND. lgauge(3,i) < 0 ) THEN
              i1=shell_ga(i)
              IF(i1 > 0) THEN
                 ptot=pp(i1)+ps(i1)+pr(i1)+rho*hh(i1)*abs(gravity)
                 gauge(30,i)=max(ptot,pmin)
              ENDIF
           ENDIF
        ENDDO
      ENDIF
 
#endif
      RETURN
      END