hm_read_sphio.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/.
!||====================================================================
!||    hm_read_sphio          ../starter/source/loads/sph/hm_read_sphio.F
!||--- called by ------------------------------------------------------
!||    lectur                 ../starter/source/starter/lectur.F
!||--- calls      -----------------------------------------------------
!||    ancmsg                 ../starter/source/output/message/message.f
!||    get_u_geo              ../starter/source/user_interface/uaccess.F
!||    hm_get_floatv          ../starter/source/devtools/hm_reader/hm_get_floatv.f
!||    hm_get_floatv_dim      ../starter/source/devtools/hm_reader/hm_get_floatv_dim.F
!||    hm_get_intv            ../starter/source/devtools/hm_reader/hm_get_intv.F
!||    hm_option_read_key     ../starter/source/devtools/hm_reader/hm_option_read_key.F
!||    hm_option_start        ../starter/source/devtools/hm_reader/hm_option_start.F
!||    subrotpoint            ../starter/source/model/submodel/subrot.F
!||    udouble                ../starter/source/system/sysfus.F
!||    usr2sys                ../starter/source/system/sysfus.F
!||--- uses       -----------------------------------------------------
!||    hm_option_read_mod     ../starter/share/modules1/hm_option_read_mod.f
!||    message_mod            ../starter/share/message_module/message_mod.f
!||    r2r_mod                ../starter/share/modules1/r2r_mod.f
!||    submodel_mod           ../starter/share/modules1/submodel_mod.F
!||====================================================================
      SUBROUTINE hm_read_sphio(ISPHIO   ,VSPHIO   ,IPART    ,IGRSURF  ,
     .                         NOD2SP   ,IPARTSP  ,ITAB     ,X        ,
     .                         MFI      ,LWASPIO  ,ITABM1   ,UNITAB   ,
     .                         LSUBMODEL,RTRANS   ,NRTRANS )
C-----------------------------------------------
C   M o d u l e s
C-----------------------------------------------
      USE message_mod
      USE r2r_mod
      USE groupdef_mod
      USE hm_option_read_mod
      USE submodel_mod    
      USE unitab_mod
      USE names_and_titles_mod , ONLY : nchartitle, ncharkey
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      "com04_c.inc"
#include      "sphcom.inc"
#include      "units_c.inc"
#include      "scr17_c.inc"
#include      "r2r_c.inc"
C-----------------------------------------------
C   D u m m y   A r g u m e n t s
C-----------------------------------------------
      INTEGER ISPHIO(NISPHIO,*), IPART(LIPART1,*),
     .        NOD2SP(*),IPARTSP(*),
     .        ITAB(*),MFI,LWASPIO,ITABM1(*)
      INTEGER, INTENT(IN) :: NRTRANS
      my_real
     .   VSPHIO(*),X(3,*)
      TYPE (SURF_)   , DIMENSION(NSURF)   :: IGRSURF
      TYPE (UNIT_TYPE_),INTENT(IN) ::UNITAB 
      TYPE(SUBMODEL_DATA) LSUBMODEL(*)
      my_real, DIMENSION(NTRANSF,NRTRANS),INTENT(IN) ::  rtrans
C-----------------------------------------------
C   L o c a l   V a r i a b l e s
C-----------------------------------------------
      INTEGER I,  N, ID, J, IDS, IDPRT, IPRT, INOD, ICEL,
     .  IPROP, ISU, NSEG,
     .  MFITMP, IVAD, LVAD, ITYPE, IDSURF, IFTEMP,
     .  ifvits,ifdens, ifpres, ifener,skip,
     .  iad,in1,in2,in3,in4, iun,nbox,nboy,nboz,nband,sub_id
      CHARACTER MESS*40
      CHARACTER(LEN=NCHARTITLE) :: TITR
      CHARACTER(LEN=NCHARKEY) :: KEY
      my_real
     .   bid,rhoin,pin,ein,dist,
     .   x1,x2,x3,x4,y1,y2,y3,y4,z1,z2,z3,z4,
     .   x12,y12,z12,x13,y13,z13,nn,nx,ny,nz,
     .   dbucs,xbmin,ybmin,zbmin,xbmax,ybmax,zbmax,
     .   pinfini,carl,xx(9),fcut,
     .   rhoin_unit,pin_unit,ein_unit,xl(3)
      LOGICAL IS_AVAILABLE,FOUND
C-----------------------------------------------
      DATA MESS/'SPH INLET/OUTLET DEFINITION             '/
      DATA iun/1/
C-----------------------------------------------
C   E x t e r n a l   F u n c t i o n s
C-----------------------------------------------
      INTEGER USR2SYS   
      my_real
     .         GET_U_GEO
      EXTERNAL get_u_geo
C-----------------------------------------------
C                                                      TYPE=1    TYPE=2   TYPE=3   TYPE=4
C   ISPHIO(1,N) = CONDITION TYPE                          x         x        x        x
C      1: INLET  2:GENERAL OUTLET  3:CONTINUE OUTLET
C   ISPHIO(2,N) = NODES GROUP RELATED TO THE CONDITION    x         x        x        x
C   ISPHIO(3,N) = SURFACE                                 x         x        x
C   ISPHIO(4,N) =                                         x         x        x
C                 IVAD : INDEX OF REAL CORRESPONDING ARRAY INTO VSPHIO
C   ISPHIO(5,N) = DENSITY FUNCTION NUMBER                 x         x                   
C   ISPHIO(6,N) = PRESSURE FUNCTION NUMBER                          x        x          
C   ISPHIO(7,N) = ENERGY FUNCTION NUMBER                  x         x                   
C   ISPHIO(8,N) = NORMAL VELOCITY FUNCTION NUMBER         x
C   ISPHIO(9,N) = CONDITION IDENTIFIER                    x         x        x        x
C   ISPHIO(10,N) = NB SEG IO                              x         x        x
C   ISPHIO(11,N) = IAD for // outlet                      x         x        x        x
C   ISPHIO(12,N) = type of input                                    x        x        x
C   ISPHIO(13,N) = id of nod1                                       x        x        x
C   ISPHIO(14,N) = id of nod2                                       x        x        x
C   ISPHIO(15,N) = id of nod3                                       x        x        x
C-----------------------------------------------
C                                                      TYPE=1    TYPE=2   TYPE=3   TYPE=4
C   VSPHIO(IVAD  ) = RHOIN                                 x         
C   VSPHIO(IVAD+1) = PIN                                            x        x        x
C   VSPHIO(IVAD+2) = EIN                                   x         
C   VSPHIO(IVAD+3) = WITHIN/WITHOUT DISTANCE               x        x        x        x      
C
Case TYPE=1
C   VSPHIO(IVAD+4+2*J,J=0,NSEG-1) =                        x
C          MASS FLOW THROUGH SEGMENT J SINCE LAST ENTRY
C   VSPHIO(IVAD+5+2*J,J=0,NSEG-1) = SEGMENTS WIDTH         x                       
Case TYPE=2,3,4 control section + new outlet / silent boundary (if surface defined by coordinates)                          
C   VSPHIO(IVAD+4->6) = X1                                          x        x        x                           
C   VSPHIO(IVAD+7->9) = X2                                          x        x        x                          
C   VSPHIO(IVAD+10->12) = X3                                        x        x        x
Case TYPE=2,3,4 control section + new outlet / silent boundary
C   VSPHIO(IVAD+13) = Total mass that cross section                 x        x        x
C   VSPHIO(IVAD+14) = DT12                                          x        x        x
C   VSPHIO(IVAD+15) = Fcut Butterworth filter                       x        x        x
C   VSPHIO(IVAD+16) = Filtered mass that cross section              x        x        x
C   VSPHIO(IVAD+17->21) = Butterworth data                          x        x        x                            
C-----------------------------------------------
      lwaspio = 0
      nseg_io = 0
      is_available = .false.
C
      mfitmp = mfi
      ivad = 1
C
      ! Start reading the option
      CALL hm_option_start('/SPH/INOUT') 
C
      ! Loop over /SPH/INOUT
      i = 0
      DO n = 1,nsphio
        skip = 0
CC --------- Multidomatic-> We do not know the inlets non-tages ------------------------------------
        IF (nsubdom > 0) THEN
          IF (tagsphio(n) == 0) skip = 1
        ENDIF
C----------------------------------------------------------------------
        IF (skip == 0) THEN
C
          ! Title and ID
          titr = ''   
          CALL hm_option_read_key(lsubmodel, 
     .                            option_id   = id, 
     .                            option_titr = titr, 
     .                          submodel_id   = sub_id)
C        
          i = i+1
          isphio(nisphio,i)=id
C
          ! Read card 1
          CALL hm_get_intv('Itype'   ,itype  ,is_available,lsubmodel)
          CALL hm_get_intv('pid'     ,idprt  ,is_available,lsubmodel)
          CALL hm_get_intv('SURF_ID' ,idsurf ,is_available,lsubmodel)
          CALL hm_get_floatv('DIST'  ,dist   ,is_available,lsubmodel, unitab)
          
          CALL hm_get_intv('node_ID1',in1 ,is_available,lsubmodel)
          CALL hm_get_intv('node_ID2',in2 ,is_available,lsubmodel)
          CALL hm_get_intv('node_ID3',in3 ,is_available,lsubmodel)
          CALL hm_get_floatv('Fcut'  ,fcut   ,is_available,lsubmodel, unitab)      
          isphio(1,i) =itype
C
          ids = 0
          found = .false.
          DO j = 1,npart
            IF (ipart(4,j) == idprt) THEN
              ids = j
              found = .true.
              EXIT
            ENDIF
          ENDDO
          IF (.NOT.found) THEN
            CALL ancmsg(msgid=437,
     .                  msgtype=msgerror,
     .                  anmode=aninfo,
     .                  i1=id,
     .                  c1=titr,
     .                  i2=idprt)
          ENDIF
          isphio(2,i) = ids
          isphio(4,i) = ivad
C
          IF (idsurf > 0) THEN
C         input by surfid
            ids=0
            found = .false.
            DO j=1,nsurf
              IF (igrsurf(j)%ID == idsurf) THEN
                ids=j
                found = .true.
                EXIT
              ENDIF
            ENDDO
            IF (.NOT.found) THEN
              CALL ancmsg(msgid=438,
     .                msgtype=msgerror,
     .                anmode=aninfo,
     .                i1=id,
     .                c1=titr,
     .                i2=idsurf)
            ENDIF
            isphio(3,i) =ids
          ELSEIF ((in1 == 0).AND.(in2 == 0).AND.(in3 == 0)) THEN
C         input by coordinates ->  ISPHIO(12,I) = 1
            isphio(12,i) =1
            CALL hm_get_floatv('XM'  ,xx(1)   ,is_available,lsubmodel, unitab)   
            CALL hm_get_floatv('YM'  ,xx(2)   ,is_available,lsubmodel, unitab)  
            CALL hm_get_floatv('ZM'  ,xx(3)   ,is_available,lsubmodel, unitab)  
            CALL hm_get_floatv('xm1' ,XX(4)   ,IS_AVAILABLE,LSUBMODEL, UNITAB)   
            CALL HM_GET_FLOATV('ym1' ,XX(5)   ,IS_AVAILABLE,LSUBMODEL, UNITAB)  
            CALL HM_GET_FLOATV('zm1' ,XX(6)   ,IS_AVAILABLE,LSUBMODEL, UNITAB)              
            CALL HM_GET_FLOATV('xm2' ,XX(7)   ,IS_AVAILABLE,LSUBMODEL, UNITAB)   
            CALL HM_GET_FLOATV('ym2' ,XX(8)   ,IS_AVAILABLE,LSUBMODEL, UNITAB)  
            CALL HM_GET_FLOATV('zm2' ,XX(9)   ,IS_AVAILABLE,LSUBMODEL, UNITAB)              
            IF(SUB_ID /= 0) THEN
              XL(1:3) = XX(1:3)
              CALL SUBROTPOINT(XL(1),XL(2),XL(3),RTRANS,SUB_ID,LSUBMODEL)
              XX(1:3) = XL(1:3) 
              XL(1:3) = XX(4:6)
              CALL SUBROTPOINT(XL(1),XL(2),XL(3),RTRANS,SUB_ID,LSUBMODEL)
              XX(4:6) = XL(1:3) 
              XL(1:3) = XX(7:9)
              CALL SUBROTPOINT(XL(1),XL(2),XL(3),RTRANS,SUB_ID,LSUBMODEL)
              XX(7:9) = XL(1:3) 
            END IF
            DO J = 1,9
              VSPHIO(IVAD+3+J) = XX(J)
            ENDDO 
          ELSE
C         input by node ID ->  ISPHIO(12,I) = 2
            ISPHIO(12,I) = 2
            ISPHIO(13,I) = USR2SYS(IN1,ITABM1,MESS,ID)
            ISPHIO(14,I) = USR2SYS(IN2,ITABM1,MESS,ID)
            ISPHIO(15,I) = USR2SYS(IN3,ITABM1,MESS,ID)  
          ENDIF
C
          IF (ITYPE == 1) THEN
            NSEG=IGRSURF(IDS)%NSEG
            LVAD=4+2*NSEG
          ELSE
            LVAD=22
          ENDIF
C
          IF (ITYPE == 1)THEN
            CALL HM_GET_INTV('fun_a1'  ,IFDENS ,IS_AVAILABLE,LSUBMODEL)         
            CALL HM_GET_FLOATV('r0k0'  ,RHOIN  ,IS_AVAILABLE,LSUBMODEL,UNITAB) 
.AND.            IF ((RHOIN == ZERO)(IFDENS>0)) THEN
              CALL HM_GET_FLOATV_DIM('r0k0' ,RHOIN_UNIT ,IS_AVAILABLE, LSUBMODEL, UNITAB)
              RHOIN = ONE * RHOIN_UNIT
            ENDIF
            CALL HM_GET_INTV('fun_a6'  ,IFENER ,IS_AVAILABLE,LSUBMODEL)         
            CALL HM_GET_FLOATV('mat_e0',EIN    ,IS_AVAILABLE,LSUBMODEL,UNITAB)  
.AND.            IF ((EIN == ZERO)(IFENER>0)) THEN
              CALL HM_GET_FLOATV_DIM('mat_e0'  ,EIN_UNIT ,IS_AVAILABLE, LSUBMODEL, UNITAB)
              EIN = ONE * EIN_UNIT
            ENDIF
            CALL HM_GET_INTV('fun_a3'  ,IFVITS ,IS_AVAILABLE,LSUBMODEL)    
          ELSEIF (ITYPE == 2) THEN
            CALL HM_GET_INTV('fun_a2'   ,IFPRES ,IS_AVAILABLE,LSUBMODEL)         
            CALL HM_GET_FLOATV('mat_p0' ,PIN    ,IS_AVAILABLE,LSUBMODEL,UNITAB) 
.AND.            IF ((PIN == ZERO)(IFPRES>0)) THEN
              CALL HM_GET_FLOATV_DIM('mat_p0'  ,PIN_UNIT ,IS_AVAILABLE, LSUBMODEL, UNITAB)
              PIN = ONE * PIN_UNIT
            ENDIF
          ELSEIF (ITYPE == 3) THEN
            CALL HM_GET_INTV('fun_a4'   ,IFPRES ,IS_AVAILABLE,LSUBMODEL)         
            CALL HM_GET_FLOATV('mat_pscale',PIN ,IS_AVAILABLE,LSUBMODEL, UNITAB) 
.AND.            IF ((PIN == ZERO)(IFPRES>0)) THEN
              CALL HM_GET_FLOATV_DIM('mat_pscale' ,PIN_UNIT ,IS_AVAILABLE, LSUBMODEL, UNITAB)
              PIN = ONE * PIN_UNIT
            ENDIF
            CALL HM_GET_FLOATV('lc'     ,CARL   ,IS_AVAILABLE,LSUBMODEL, UNITAB)     
          ENDIF
C
          IF (ITYPE == 1) THEN
           ISPHIO(5,I) = IFDENS
           ISPHIO(7,I) = IFENER
           ISPHIO(8,I) = IFVITS
          ELSEIF (ITYPE == 2) THEN
           ISPHIO(6,I) = IFPRES
          ELSEIF (ITYPE == 3) THEN
           ISPHIO(6,I) = IFPRES
          ENDIF
!
          IF (ITYPE == 1) THEN
.AND.            IF ((IFDENS /= 0)(RHOIN == ZERO)) RHOIN = ONE
.AND.            IF ((IFENER /= 0)(EIN == ZERO))   EIN   = ONE
            VSPHIO(IVAD  )  = RHOIN
            VSPHIO(IVAD+2)  = EIN
            VSPHIO(IVAD+3)  = DIST
          ELSEIF (ITYPE == 2) THEN
            VSPHIO(IVAD+1)  = PIN
            VSPHIO(IVAD+3)  = DIST
            VSPHIO(IVAD+15) = FCUT
          ELSEIF (ITYPE == 3) THEN
            VSPHIO(IVAD+1)  = PIN
            VSPHIO(IVAD+2)  = CARL
            VSPHIO(IVAD+3)  = DIST
            VSPHIO(IVAD+15) = FCUT
          ELSEIF (ITYPE == 4) THEN
            VSPHIO(IVAD+15) = FCUT
          ENDIF
          MFITMP = MFITMP+LVAD
          IVAD   = IVAD+LVAD
C----------------------------------------------------------------------
        ENDIF
      ENDDO
      MFI = MFITMP
C-------------------------------------
C search for duplicate ids
C-------------------------------------
      CALL UDOUBLE(ISPHIO(NISPHIO,1),NISPHIO,NSPHIO,MESS,0,BID)
C-------------------------------------
C inlets, pre-computes surface of segments.
C-------------------------------------
      DO N=1,NSPHIO
       IF(ISPHIO(1,N)==1)THEN
        IVAD=ISPHIO(4,N)+5
        ISU =ISPHIO(3,N)
        NSEG=IGRSURF(ISU)%NSEG
        DO J=1,NSEG
          IN1=IGRSURF(ISU)%NODES(J,1)
          IN2=IGRSURF(ISU)%NODES(J,2)
          IN3=IGRSURF(ISU)%NODES(J,3)
          IN4=IGRSURF(ISU)%NODES(J,4)
          X1=X(1,IN1)
          Y1=X(2,IN1)
          Z1=X(3,IN1)
          X2=X(1,IN2)
          Y2=X(2,IN2)
          Z2=X(3,IN2)
          X3=X(1,IN3)
          Y3=X(2,IN3)
          Z3=X(3,IN3)
          X12=X1-X2
          Y12=Y1-Y2
          Z12=Z1-Z2
          X13=X3-X2
          Y13=Y3-Y2
          Z13=Z3-Z2
          NX=Y12*Z13-Z12*Y13
          NY=Z12*X13-X12*Z13
          NZ=X12*Y13-Y12*X13
          NN =SQRT(NX*NX+NY*NY+NZ*NZ)
          IF(IN4/=IN3)THEN
           X4=X(1,IN4)
           Y4=X(2,IN4)
           Z4=X(3,IN4)
           X12=X1-X4
           Y12=Y1-Y4
           Z12=Z1-Z4
           X13=X3-X4
           Y13=Y3-Y4
           Z13=Z3-Z4
           NX=Y12*Z13-Z12*Y13
           NY=Z12*X13-X12*Z13
           NZ=X12*Y13-Y12*X13
           NN =NN+SQRT(NX*NX+NY*NY+NZ*NZ)
          ENDIF
          VSPHIO(IVAD)=HALF*NN
          IVAD=IVAD+2
        ENDDO
       ENDIF
      ENDDO
C-------------------------------------
C inlets, outlets : necessary work space (cf engine).
C-------------------------------------
      DO N=1,NSPHIO
        IF(ISPHIO(12,N)==0)THEN
        IVAD=ISPHIO(4,N)
        ISU =ISPHIO(3,N)
        NSEG=IGRSURF(ISU)%NSEG
        NSEG_IO = NSEG_IO + NSEG
        DIST =VSPHIO(IVAD+3)
        DBUCS=DIST
        XBMIN =1.E+20
        YBMIN =1.E+20
        ZBMIN =1.E+20
        XBMAX=-1.E+20
        YBMAX=-1.E+20
        ZBMAX=-1.E+20
        DO J=1,NSEG
          IN1=IGRSURF(ISU)%NODES(J,1)
          IN2=IGRSURF(ISU)%NODES(J,2)
          IN3=IGRSURF(ISU)%NODES(J,3)
          IN4=IGRSURF(ISU)%NODES(J,4)
          X1=X(1,IN1)
          Y1=X(2,IN1)
          Z1=X(3,IN1)
          X2=X(1,IN2)
          Y2=X(2,IN2)
          Z2=X(3,IN2)
          X3=X(1,IN3)
          Y3=X(2,IN3)
          Z3=X(3,IN3)
          XBMIN=MIN(XBMIN,X1)
          YBMIN=MIN(YBMIN,Y1)
          ZBMIN=MIN(ZBMIN,Z1)
          XBMAX=MAX(XBMAX,X1)
          YBMAX=MAX(YBMAX,Y1)
          ZBMAX=MAX(ZBMAX,Z1)
          XBMIN=MIN(XBMIN,X2)
          YBMIN=MIN(YBMIN,Y2)
          ZBMIN=MIN(ZBMIN,Z2)
          XBMAX=MAX(XBMAX,X2)
          YBMAX=MAX(YBMAX,Y2)
          ZBMAX=MAX(ZBMAX,Z2)
          XBMIN=MIN(XBMIN,X3)
          YBMIN=MIN(YBMIN,Y3)
          ZBMIN=MIN(ZBMIN,Z3)
          XBMAX=MAX(XBMAX,X3)
          YBMAX=MAX(YBMAX,Y3)
          ZBMAX=MAX(ZBMAX,Z3)
          DBUCS=MAX(DBUCS,ABS(X1-X2))
          DBUCS=MAX(DBUCS,ABS(Y1-Y2))
          DBUCS=MAX(DBUCS,ABS(Z1-Z2))
          DBUCS=MAX(DBUCS,ABS(X2-X3))
          DBUCS=MAX(DBUCS,ABS(Y2-Y3))
          DBUCS=MAX(DBUCS,ABS(Z2-Z3))
          DBUCS=MAX(DBUCS,ABS(X3-X1))
          DBUCS=MAX(DBUCS,ABS(Y3-Y1))
          DBUCS=MAX(DBUCS,ABS(Z3-Z1))
          IN4=IGRSURF(ISU)%NODES(J,4)
          IF(IN4/=IN3)THEN
           X4=X(1,IN4)
           Y4=X(2,IN4)
           Z4=X(3,IN4)
           XBMIN=MIN(XBMIN,X4)
           YBMIN=MIN(YBMIN,Y4)
           ZBMIN=MIN(ZBMIN,Z4)
           XBMAX=MAX(XBMAX,X4)
           YBMAX=MAX(YBMAX,Y4)
           ZBMAX=MAX(ZBMAX,Z4)
           DBUCS=MAX(DBUCS,ABS(X1-X4))
           DBUCS=MAX(DBUCS,ABS(Y1-Y4))
           DBUCS=MAX(DBUCS,ABS(Z1-Z4))
           DBUCS=MAX(DBUCS,ABS(X2-X4))
           DBUCS=MAX(DBUCS,ABS(Y2-Y4))
           DBUCS=MAX(DBUCS,ABS(Z2-Z4))
           DBUCS=MAX(DBUCS,ABS(X3-X4))
           DBUCS=MAX(DBUCS,ABS(Y3-Y4))
           DBUCS=MAX(DBUCS,ABS(Z3-Z4))
          ENDIF
        ENDDO
        XBMIN=XBMIN-DIST
        YBMIN=YBMIN-DIST
        ZBMIN=ZBMIN-DIST
        XBMAX=XBMAX+DIST
        YBMAX=YBMAX+DIST
        ZBMAX=ZBMAX+DIST
        NBOX =MAX(IUN,INT((XBMAX-XBMIN)/DBUCS))
        NBOY =MAX(IUN,INT((YBMAX-YBMIN)/DBUCS))
        NBOZ =MAX(IUN,INT((ZBMAX-ZBMIN)/DBUCS))
        NBAND=MAX(NBOX,NBOY,NBOZ)+1
        ELSE
          NBAND = 2
        ENDIF
        LWASPIO=MAX(LWASPIO,15*NUMSPH+6*(NBAND+1)+12*NSEG)
      ENDDO
      LWASPIO=MAX(LWASPIO,3*NSPHIO)
C-------------------------------------
C print
C-------------------------------------
      WRITE(IOUT,1000)
      DO N=1,NSPHIO
         IVAD=ISPHIO(4,N)
         IF(ISPHIO(1,N)==1)THEN
           WRITE(IOUT,1100) ISPHIO(NISPHIO,N),
     .     IPART(4,ISPHIO(2,N)),IGRSURF(ISPHIO(3,N))%ID,
     .     ISPHIO(5,N),VSPHIO(IVAD),ISPHIO(7,N),
     .     VSPHIO(IVAD+2),ISPHIO(8,N),VSPHIO(IVAD+3)
         ELSE
           IF(ISPHIO(1,N)==2)THEN
             IF (ISPHIO(12,N)==0) THEN
               WRITE(IOUT,1200) ISPHIO(NISPHIO,N),
     .         IPART(4,ISPHIO(2,N)),IGRSURF(ISPHIO(3,N))%ID,
     .         ISPHIO(6,N),VSPHIO(IVAD+1),VSPHIO(IVAD+3)
             ELSE
               WRITE(IOUT,1400) ISPHIO(NISPHIO,N),
     .         IPART(4,ISPHIO(2,N)),ISPHIO(6,N),VSPHIO(IVAD+1),VSPHIO(IVAD+3)
             ENDIF
           ELSEIF(ISPHIO(1,N)==3)THEN
             IF (ISPHIO(12,N)==0) THEN
               WRITE(IOUT,1300) ISPHIO(NISPHIO,N),
     .         IPART(4,ISPHIO(2,N)),IGRSURF(ISPHIO(3,N))%ID,VSPHIO(IVAD+3),
     .         ISPHIO(6,N),VSPHIO(IVAD+1),VSPHIO(IVAD+2)
             ELSEIF (ISPHIO(12,N)==1) THEN
               WRITE(IOUT,1500) ISPHIO(NISPHIO,N),IPART(4,ISPHIO(2,N)),
     .         ISPHIO(6,N),VSPHIO(IVAD+1),VSPHIO(IVAD+2)
             ENDIF
           ELSEIF(ISPHIO(1,N)==4)THEN
             IF (ISPHIO(12,N)==0) THEN
               WRITE(IOUT,1600) ISPHIO(NISPHIO,N),IPART(4,ISPHIO(2,N)),
     .         IGRSURF(ISPHIO(3,N))%ID
             ELSE
               WRITE(IOUT,1700) ISPHIO(NISPHIO,N),IPART(4,ISPHIO(2,N))
             ENDIF
           ENDIF
           IF (ISPHIO(12,N)==1) THEN
C            input by coordinates
             IVAD=ISPHIO(4,N)
             WRITE(IOUT,2100) VSPHIO(IVAD+4),VSPHIO(IVAD+5),VSPHIO(IVAD+6),
     .       VSPHIO(IVAD+7),VSPHIO(IVAD+8),VSPHIO(IVAD+9),VSPHIO(IVAD+10),
     .       VSPHIO(IVAD+11),VSPHIO(IVAD+12) 
           ELSEIF (ISPHIO(12,N)==2) THEN
C            input by node ID
             WRITE(IOUT,2200) ITAB(ISPHIO(13,N)),ITAB(ISPHIO(14,N)),ITAB(ISPHIO(15,N))
           ENDIF
           IF (VSPHIO(IVAD+15)>EM20) WRITE(IOUT,2300) VSPHIO(IVAD+15)
         ENDIF
      ENDDO
C-------------------------------------
C
      RETURN
C
 1000 FORMAT(//
     .'     sph inlet/outlet conditions       '/
     .'     ---------------------------       ')
 1100 FORMAT(/5X ,'sph inlet/outlet condition id           ',I10
     .       /5X ,'TYPE inlet                              ',
     .       /10X,'part related to condition               ',I10
     .       /10X,'inlet surface                           ',I10
     .       /10X,'density FUNCTION                        ',I10
     .       /10X,'scale factor on density function        ',1PG20.13
     .       /10X,'energy function                         ',I10
     .       /10X,'scale factor on energy function         ',1PG20.13
     .       /10X,'normal velocity function                ',I10
     .       /10X,'within distance for particles setting   ',1PG20.13)
 1200 FORMAT(/5X ,'sph inlet/outlet condition id           ',I10
     .       /5X ,'type general outlet                     ',
     .       /10X,'part related to condition               ',I10
     .       /10X,'outlet surface                          ',I10
     .       /10X,'pressure function                       ',I10
     .       /10X,'scale factor on pressure function       ',1PG20.13
     .       /10X,'without distance for particles setting  ',1PG20.13)
 1300 FORMAT(/5X ,'sph inlet/outlet condition id           ',I10
     .       /5X ,'type silent boundary                    ',
     .       /10X,'part related to condition               ',I10
     .       /10X,'outlet surface                          ',I10
     .       /10X,'without distance for particles setting  ',1PG20.13
     .       /10X,'pressure function                       ',I10
     .       /10X,'scale factor on pressure function       ',1PG20.13
     .       /10X,'characteristic length                   ',1PG20.13)
 1400 FORMAT(/5X ,'sph inlet/outlet condition id           ',I10
     .       /5X ,'type general outlet                     ',
     .       /10X,'part related to condition               ',I10
     .       /10X,'pressure function                       ',I10
     .       /10X,'scale factor on pressure function       ',1PG20.13
     .       /10X,'without distance for particles setting  ',1PG20.13)
 1500 FORMAT(/5X ,'sph inlet/outlet condition id           ',I10
     .       /5X ,'type silent boundary                    ',
     .       /10X,'part related to condition               ',I10
     .       /10X,'pressure function                       ',I10
     .       /10X,'scale factor on pressure function       ',1PG20.13
     .       /10X,'without distance for particles setting  ',1PG20.13
     .       /10X,'characteristic length                   ',1PG20.13)
 1600 FORMAT(/5X ,'sph inlet/outlet condition id           ',I10
     .       /5X ,'type sph control section                ',
     .       /10X,'part related to control section         ',I10
     .       /10X,'surface                                 ',I10)
 1700 FORMAT(/5X ,'sph inlet/outlet condition id           ',I10
     .       /5X ,'type sph control section                ',
     .       /10X,'part related to control section         ',I10)
C
 2100 FORMAT(10X,'surface defined by coordinates          ',
     .       /10X,' --> coordinates of node1               ',1PG20.13,1PG20.13,1PG20.13,
     .       /10X,' --> coordinates of node2               ',1PG20.13,1PG20.13,1PG20.13,
     .       /10X,' --> coordinates of node3               ',1PG20.13,1PG20.13,1PG20.13)
C
 2200 FORMAT(10X,'surface defined by nodes                ',I10,I10,I10)
C
 2300 FORMAT(10X,'4-pole butterworth corner frequency     ',1PG20.13)
C
      RETURN
      END