c

c       filename  gxpist.htm              070405
C**** SUBROUTINE GXPIST is called from GREX3, group 19, section 1,
C     and is activated only when W1AD is set to ZMOVE and when AZW1
C     is set to greater than zero in SATELLITE.
C     Click here for the
C     the PHOENICS Encyclopaedia description.
C
C.... The grid consists of two parts. The first part covers the zone
C     from the top of the cylinder ( ie. z=0.0) to the top of the
C     piston; this part contracts and expands as the piston moves
C     up and down. The second part of the grid covers the bowl of
C     the cylinder, & hence moves with the cylinder but does not
C     expand or contract. The coding which follows assumes that the
C     satellite-set grid distribution is the one appropriate for
C     "bottom dead centre ( ie. TFIRST=0.0 ), even when a restart
C     is in question, in which case TFIRST is not zero.
C
C.... The satellite-set parameters which enter here are:
C     AZW1, the rotation rate in radians/sec;
C     BZW1, the crank radius; and
C     CZW1, the ratio of the connecting-rod length to the crank radius.
C     IZW1, the iz location of the last slab in part one of an n-part
C           moving grid.
C.... The library cases 327 and 330-331 make use of it.
C
      SUBROUTINE GXPIST(ISTEP,NZ,TIM,DT,ISWEEP,LSWEEP,NPRINT,NTPRIN)
      INCLUDE 'farray'
      INCLUDE 'grdloc'
      INCLUDE 'satgrd'
      COMMON /LDATA/ LDAT(84)
      LOGICAL LDAT,NULLPR
      DIMENSION ZWO(2),ZWN(2),WAV(2),IZW(2)
      SAVE IZW,DZBOWL,ZTDC,ZWO,ZWN,WAV
      COMMON /NAMFN/NAMFUN,NAMSUB
      CHARACTER*6 NAMFUN,NAMSUB
      EQUIVALENCE (NULLPR,LDAT(32))
C
      NAMSUB = 'GXPIST'
c.... precaution against incorrect nz setting
      IF(NZ.LT.IZW1) RETURN
      IF(GRNDNO(1,AZW1)) GO TO 100
      IF(ISTEP.EQ.1) THEN
c.... initialization
        IF(BZW1.LE.0.0.OR.CZW1.LE.0.0) THEN
          CALL WRIT40('unsuitable settings for piston movement')
          CALL SET_ERR(440,
     *        'Error. Unsuitable settings for piston movement.',1)
C          STOP
        ENDIF
c.... the iz values of the ends of the first and second parts of the grid
        IZW(1) = IZW1
        IZW(2) = NZ
c
C.... ZWNZ is the EARTH index for distances of the HIGH faces from
C     the z=0.0 plane.
c
C.... Function VARZ(INDVAR,IZZ) yields the F-array value, which related
C     to z-wise variables, corresponding to INDVAR.
c.... the depth of the bowl
        DZBOWL = VARZ(ZWNZ,NZ) - VARZ(ZWNZ,IZW1)
c
c.... the z value at "top-dead-centre", where part 1 of the grid has minimum
c     depth
        ZTDC = VARZ(ZWNZ,IZW1) - 2.0*BZW1
c
c.... the crank angle at the beginning of the first time step
        ANGLE = AZW1 * (TIM-DT)
        TERM = (SIN(ANGLE)/CZW1)**2
c
c.... the piston face position the
        ZWO(1) = ZTDC + BZW1 * (1.+COS(ANGLE)+CZW1*(1.-SQRT(1.- TERM)))
c
c.... the bottom of bowl position then
        ZWO(2) = ZWO(1) + DZBOWL
      ENDIF
C.... conditions at the end of the current time step
      ANGLE = AZW1*TIM
      TERM = (SIN(ANGLE)/CZW1)**2
c
c.... new positions of piston face and bowl bottom
      ZWN(1) = ZTDC + BZW1* (1.+COS(ANGLE)+CZW1*(1.-SQRT(1.- TERM)))
      ZWN(2) = ZWN(1) + DZBOWL
c
c.... corresponding velocities
      WAV(1) = (ZWN(1)-ZWO(1))/DT
      WAV(2) = WAV(1)
      GO TO 10
  100 CONTINUE
c.................................. fixed piston velocity = bzw1
      IF(ISTEP.EQ.1) THEN
        IZW(1) = IZW1
        IZW(2) = NZ
        ZWN(2) = VARZ(ZWNZ,NZ)
        ZWN(1) = VARZ(ZWNZ,IZW1)
        WAV(1) = BZW1
        WAV(2) = BZW1
      ENDIF
      ZWO(2) = ZWN(2)
      ZWO(1) = ZWN(1)
      ZWN(2) = ZWN(2)  + BZW1*DT
      ZWN(1) = ZWN(1)  + BZW1*DT
 
   10 CALL ZMOVE(IZW,WAV,ZWN,ZWO,2)
c.... print-out
      IF(ISWEEP.EQ.LSWEEP .OR. MOD(ISWEEP,NPRINT).EQ.0) THEN
C.... ZWN is the grid location at the present time
C     ZWO is the grid location at the previous time
C     WAV is the velocity at which the grid moves
        IF(.NOT.NULLPR) THEN
          IF(MOD(ISTEP,NTPRIN).EQ.0) THEN
            CALL WRIT40('two-part expanding-grid parameters at:')
            CALL WRIT1R('time    ',tim)
            CALL WRIT2R('new zpis',zwn(1),'new zbwl',zwn(2))
            CALL WRIT2R('old zbwl',zwo(1),'old zbwl',zwo(2))
            CALL WRIT2R('vel pist',wav(1),'vel bowl',wav(2))
          ENDIF
        ENDIF
      ENDIF
      NAMSUB = 'gxpist'
      END
c