C TS 07/2000 Possible Cseal<0
C calcsolnm.f> hseal.f : 04/91 for moment coefficients
C
C
C
C       ZEROTH ORDER SOLUTION: FLOW, FORCES & TORQUE
C                for an Annular pressure seal
C
C
c  ####     ##    #        ####    ####    ####   #       #    #  #    #
c #    #   #  #   #       #    #  #       #    #  #       ##   #  ##  ##
c #       #    #  #       #        ####   #    #  #       # #  #  # ## #
c #       ######  #       #            #  #    #  #       #  # #  #    #
c #    #  #    #  #       #    #  #    #  #    #  #       #   ##  #    #
c  ####   #    #  ######   ####    ####    ####   ######  #    #  #    #


C *****************************************************************************
C **                                                                         **
C **  SETBC:  Imposses known boundary conditions on closure of domain.       **
C **                                                                         **
C *****************************************************************************

      SUBROUTINE SETBC

      IMPLICIT NONE

      INCLUDE 'params.f' 

C ----------------------------------------------------------------------------
C --  COMMON variable declarations                                         --
C ----------------------------------------------------------------------------

      COMMON /UVARRAY/ U(MAXNXT,MAXNYI), V(MAXNXT,MAXNYI)
      COMMON /PARRAY/ P(MAXNXT,0:MAXNYI)
      COMMON /FACTORS/ REP, REY, MSPEED, SPEED,
     +                 KLOSY, SWIRL, PCAV, ALFU, BETU, ALFP
      common /param3/emuc,rhoc,rpm,Ps,Pa,Pc,Losxsi,Alpha
      COMMON /ENDSEAL/ CSEAL, CDSEAL
      COMMON /NODES/ Nx,Ny

      DOUBLE PRECISION U, V, P
      DOUBLE PRECISION Pext, REP,REY,MSPEED,SPEED,CSEAL,CDSEAL,
     +                 KLOSY, SWIRL, PCAV, ALFU, BETU, ALFP,
     +                 emuc,rhoc,rpm,Ps,Pa,Pc,Losxsi,Alpha
      INTEGER Nx, Ny, j

      
C ----------------------------------------------------------------------------
C --  SETBC code                                                            --
C ----------------------------------------------------------------------------
      IF (Ps.eq.Pa) THEN
         Pext=0.0D0
      ELSE
         Pext=1.0D0 
      END IF

      DO J=1, NX                       !
          P(J, 0)=Pext                 ! External Psupply at y=0
          P(J, 1)=Pext                 !
          U(J, 1)=SWIRL                ! inlet circumferential speed
      END DO                           !
      IF (CSEAL.eq.0.0D0) THEN
        DO J=1, NX
          P(J, NY)=0.0D0               ! Zero differential end pressure at y=L/R
        END DO
      END IF                                       
      END 


C *****************************************************************************
C **                                                                         **
C **  Subroutine Guesp                                                       **
C **                                                                         **
C **  GUESP:  Guess of pressure field for concentric operation               **
C **                                                                         **
C *****************************************************************************

      SUBROUTINE GUESP(DEVICE)

      IMPLICIT NONE

      INCLUDE 'params.f'

C ----------------------------------------------------------------------------
C --  COMMON variable declarations                                         --
C ----------------------------------------------------------------------------

      COMMON /XYVEC/ XP(MAXNXT), XU(MAXNXT), YP(MAXNYI), YV(MAXNYI)
      COMMON /UVARRAY/ U(MAXNXT,MAXNYI), V(MAXNXT,MAXNYI)
      COMMON /PARRAY/ P(MAXNXT,0:MAXNYI)
      COMMON /HFILM/ hp(maxnxt,maxnyi), hu(maxnxt,maxnyi),
     +               hv(maxnxt,maxnyi)

      COMMON /PARAM3/ emuc,rhoc,rpm,Ps,Pa,Pc,Losxsi,Alpha
      COMMON /FACTORS/ REP, REY, MSPEED, SPEED,
     +                 klosy, SWIRL, PCAV, ALFU, BETU, ALFP
      COMMON /ENDSEAL/ CSEAL, CDSEAL
      COMMON /Flow0/ Qin, Qout, Qfactor 
      common /zeropars/ L, Vin, Fy, Pedge

      COMMON /FLAGS/ INERL, INERP, ITURB, ICAV
      COMMON /NODES/ Nx, Ny
      COMMON /VERB/ SVERB, DVERB, BEEP
C    !................................................!    
      DOUBLE PRECISION XP, XU, YP, YV, hp, hu, hv,
     +                 U, V, P, Qin, QOut, Qfactor,
     +                 emuc, rhoc, rpm, Ps, Pa, Pc, losxsi, 
     +                 alpha, REP, REY, MSPEED, SPEED,
     +                 klosy, SWIRL, PCAV, ALFU, BETU, ALFP,
     +                 CSEAL, CDSEAL
      INTEGER Nx, Ny,INERL,INERP,ITURB,ICAV,SVERB,DVERB,BEEP

C ----------------------------------------------------------------------------
C --  Local variable declarations                                          --
C ----------------------------------------------------------------------------

      DOUBLE PRECISION KX, KY, UIN, Y, REAX, VIN, Schilds,
     +                 PEDGE, Hi, Vold, Fy, Fz, Send, 
     +                 FNKY, FNKX, Kxrot, Dumy, L, B, Ue, Uee, 
     +                 pratio, Ay, Emu, Rho, Pj, PI, Mflow,
     +                 Rhout, Emout, Vout, Pout

      INTEGER J, K, KP1, I, DEVICE

C ----------------------------------------------------------------------------
C --  GUESP code                                                            --
C ----------------------------------------------------------------------------
      IF (Ps.eq.Pa) THEN      ! FOR PLAIN FLOODED JOURNAL BEARING
         Uin=Mspeed
         Vin=0.0D0
         Pedge=0.0D0
         KY=FNKY(Uin,Vin,1.0D0, SPEED, REP, RHO, EMU)
         Alpha=0.5D0
         Swirl=Mspeed
         DO  I=1, Nx 
          DO J=1,Ny                   
            U(i,j)=Mspeed
            V(i,j)=Vin                 
            P(i,j)=0.0D0              
          END DO   
            P(i,0)=0.0D0                  
      END DO                        

      RETURN

      END IF
c ----------------------------------------!-----------------------------   

      Pratio=1.0D0                        != Psupply

 200  CALL LOCPROPS(RHO,EMU,PRATIO)       !
      CALL LOCPROPS(Rhout,Emout,0.0D0)
      KX=12.0D0*EMU                       ! Laminar shear forces
      KY=12.0D0*EMU                       !
      Ay=Klosy*Inerp                      !
      Hi=Hp(1,1)                          !                                          ! 
      UIN=SWIRL                           ! Circ. mean velocity
      Y=YP(Ny)                            !

c      IF(Inerl.eq.0) tHEN                 !
c         Uin=Mspeed                       !
c         Alpha=0.5                        !
c         Swirl=Mspeed                     !
c      END IF                              !

      REAX=REP*PRATIO*(RHO/EMU)/Y         ! Axial press reynolds number

      KY=.066D0*REAX**.75                 ! Use Hirs formulae for approx
      KY=KY**(1.0/1.75D0)                 ! Shear correction factor
      DUMY=12.0D0*EMU                     !
      KY=DMAX1(DUMY, KY)                  ! Mean axial velocity
      VIN=PRATIO/(KY*Y)                   !
      I=0

 201  DUMY=1.0D0                          !
      KY=FNKY(UIN, VIN,Hi, SPEED, REP, RHO, EMU)
      I=I+1                               !
      Fy=DABS(Ky/Vin)                     !
      Vold=Vin                            !
      Vin=1.D0/DSQRT(Fy*Y+RHO*(Ay+CDSEAL*RHO/Rhout))       
      PEDGE=PRATIO-AY*VIN*VIN             ! Calculate approximate value
      call locprops(Rho,Emu,Pedge)        ! of entrance pressure
      Reax=Rep*(Rho/Emu)*Vin/Y            !
      Vout=Vin*RHO/Rhout
      Pout=Cdseal*Vout*Vout*Rhout
      call locprops(Rhout,Emout,Pout)

      IF (DABS((Vold-Vin)/Vin).le.(0.001)) then
c    !....................................!
      IF (SVERB) THEN                     !
          WRITE (6, 20) PRATIO, VIN, PEDGE, KY
      END IF                              !
      IF (DEVICE.AND.(DVERB.EQ.3)) THEN   !
          WRITE (1, 20) PRATIO, VIN, PEDGE, KY
      END IF                              !

C    !....................................!

      IF (PEDGE.LE.PCAV) THEN             ! Warning if Pedge<Pcav
          WRITE (6, *)                    !
     +       'GUESSED ENTRANCE Pressure TOO LOW,',
     +       'Loss factor too large, PRG STOPS'
          CALL BEEPER                     !
          STOP
      END IF                              !
C    !....................................!
      ELSE
       goto 201
      END IF
c    !....................................! 

c...........................................................................
c Approximate soln. for concentric operation
c short seal, non-tapered solution
c...........................................................................
      PI=DACOS(-1.0D0)
      Schilds=Fy*Y/Rey                    ! Sigma coeffic. of Childs

      Kx=FNKX(Uin,Vin,Hi,Speed,Rep,Rho,Emu,Kxrot)
      L=Kx/Rey/Rho/Vin                    !
      B=(2.0D0*Alpha-1.0D0)*Inerl         !

      Ue=Mspeed*(1.0D0+B*DEXP(-L*Y))      ! exit circ. velocity
      IF (Speed.gt.(0.0D0)) THEN          !...!
          Uee=Ue/Speed                        !
      ELSE                                    !
          Uee=0.0D0                           !
      END IF                              !...!

      Qin=2.0D0*Vin*PI*Hi*Rho             ! dim. mass flow
      Mflow=Qin*Qfactor                   !


      IF (SVERB) THEN                     !
        write(6,29) Vin, Reax, Ky, Kx,Pedge,Alpha,Uee,Mflow,Schilds         
      END IF                              !
      IF (DEVICE.AND.(DVERB.EQ.3)) THEN   !
        write(1,29) Vin, Reax, Ky,Kx, Pedge, Alpha,Uee,Mflow,Schilds
      END IF                              !


      DO J=1, Ny                          ! Zero-th order guess fields
          Pj=Pedge*(1.0D0-Yp(j)/Y)+Pout*Yp(j)/Y
          Fz=DEXP(-L*Yp(j))               !
          DO  I=1, Nx                     !
            U(i,j)=Mspeed*(1.0D0+B*Fz)    !
            V(i,j)=Vin                    !
            P(i,j)=Pj                     !
          END DO                          !
      END DO                              !

      DO j=1, Nx                          ! at inlet/ exit boundary 
          P(J, NY)=Pout                   ! p=0.0
          P(J,1) = Pedge                  !
      END DO                              !
C
C.........................................!.........................

                                          
   20 FORMAT (' ', /, 4X, 'ENtrance pressure ratio, Pe/Ps:',
     +        E12.5E2, /, 4X, 'GUESSED: Vinlet=', E12.5E2, 3X,
     +        'Pedge=', E12.5E2, /, 4X, 'Shear Factor Ky=',
     +        E12.5E2)
   29 FORMAT (' ', /, 70('.'),/,4X,
     + 'Vaxial:', E12.5E2, 2X, 'Axial Rea:', E12.5E2, /, 4X,
     + 'Shear Ky:', E12.5E2, 2X, 'Ky:', E12.5E2, /, 4X,
     + 'Entrance Pressure Pe/Ps:', E12.5E2, /,4X,
     + 'Swirl/wR:', E12.5E2, 2X, 'Uexit/wR:', E12.5E2, /,
     + 'Mass Flow:', E12.5E2, ' Kg/s', /, 4X,
     + 'Friction factor a:', E12.5E2, /,
     + ' ', 70('.'),/)

      END 


C *****************************************************************************
C **                                                                         **
C **  Subroutine Sland                                                       **
C **                                                                         **
C **  SLAND:  Iterative solution of flow field on land-sills regions.        **
C **                                                                         **
C *****************************************************************************


      SUBROUTINE SLAND(DEVICE)
	use DFWIN

      IMPLICIT NONE

      INCLUDE 'params.f'

C ----------------------------------------------------------------------------
C --  COMMON variable declarations                                         --
C ----------------------------------------------------------------------------

      COMMON /DXVEC/ DXP(MAXNXT), DXU(MAXNXT), SUW(MAXNXT), SUE(MAXNXT)
      COMMON /DYVEC/ DYP(MAXNYI), DYV(MAXNYI), SVN(MAXNYI), SVS(MAXNYI)
      COMMON /HFILM/ HP(MAXNXT,MAXNYI), HU(MAXNXT,MAXNYI),
     +               HV(MAXNXT,MAXNYI) 
      COMMON /UVARRAY/ U(MAXNXT,MAXNYI), V(MAXNXT,MAXNYI)
      COMMON /PARRAY/ P(MAXNXT,0:MAXNYI)
      COMMON /RHOEMU/ RHOP(MAXNXT,0:MAXNYI), EMUP(MAXNXT,0:MAXNYI)
      COMMON /DPUV/ DU(MAXNXTP2), DV(MAXNXTP2), DVV(MAXNXTP2)
      COMMON /TDMA0/ A(MAXNXTP2), B(MAXNXTP2), C(MAXNXTP2), D(MAXNXTP2)

      COMMON /FACTORS/ REP, REY, MSPEED, SPEED,
     +                 Klosy, SWIRL, PCAV, ALFU, BETU, ALFP
      COMMON /ENDSEAL/ CSEAL, CDSEAL
      COMMON /FLOW0/ Qin, Qout, Qfactor
      COMMON /SOURCEA/ SMASS, PEPS, MMP, PMAXP
      COMMON /KVALA/ DYVK, DYPK, SVNK, SVSK

      COMMON /SOURCEB/ ITER, ITMAX, ITPMAX
      COMMON /UVEC/ JUMIN, JUMAX, JUSTART, JUSTOP
      COMMON /VVEC/ JVMIN, JVMAX, JVSTART, JVSTOP
      COMMON /PRES/ JPMIN, JPMAX, JPSTART, JPSTOP
      COMMON /FLAGS/ INERL, INERP, ITURB, Icav
      COMMON /NODES/ Nxt, Nyi
      COMMON /KVALB/ K, KM1, KP1
      COMMON /VERB/ SVERB, DVERB, iBEEP
c.........................................................................

      DOUBLE PRECISION DXP, DXU,SUW, SUE, HP, HU, HV,
     +                 DYP, DYV, SVN, SVS,
     +                 U, V, P, RHOP, EMUP,
     +                 DU, DV, DVV,A, B, C, D
      DOUBLE PRECISION QIN, QOUT, QFACTOR, CSEAL, CDSEAL,
     +                 REP, REY, MSPEED, SPEED,
     +                 klosy, SWIRL, PCAV, ALFU, BETU, ALFP,
     +                 SMASS,PEPS,MMP,PMAXP,DYVK,DYPK,SVNK,SVSK
       
      INTEGER JUMIN, JUMAX, JUSTART, JUSTOP,
     +        JVMIN, JVMAX, JVSTART, JVSTOP,
     +        JPMIN, JPMAX, JPSTART, JPSTOP,
     +        Nxt, Nyi, INERL, INERP, ITURB, ICAV,
     +        K, KM1, KP1, ITER, ITMAX, ITPMAX,
     +        SVERB, DVERB, iBEEP

C ----------------------------------------------------------------------------
C --  Local variable declarations                                          --
C ----------------------------------------------------------------------------

      DOUBLE PRECISION Sflow, Pmax, relaxp, relaxu, Sold
      INTEGER Inerldum,Dumyi,I,J,DEVICE,countp,maxcountp,iconv
 	integer(SINT) :: ret 

C ----------------------------------------------------------------------------
C --  SLAND code                                                            --
C ----------------------------------------------------------------------------
      Inerldum=Inerl    
      Sflow=Pmaxp
      relaxp=alfp
      relaxu=alfu
      iconv=0
c    !................................ Calculate dim. liquid props.
      do k=0, Nyi
         do j=1, Nxt
            call locprops(Rhop(j,k),Emup(j,k),P(J,k))
         end do
      end do
c    !...................................!............................!
 
      CALL SETLIM                        ! Sets limits for sweeps

      ITER=0                             ! Start of iterative process
      Countp=0                           !
      maxcountp=(Nxt-1)*(Nyi-1)          ! MAX internal unknown Ps
      IF (CDSEAL.gt.0.0D0) maxcountp=maxcountp+Nxt-1

  100 ITER=ITER+1                        ! Zeroes calculation vectors
C    !------------------!                !---------------------------

      DO I=1, NXT+2                        !
          DVV(I)=0.0D0                     !
          A(I)=0.0D0                       !
          B(I)=0.0D0                       !
          C(I)=0.0D0                       !
          D(I)=0.0D0                       !
          DU(I)=0.0D0                      !
          DV(I)=0.0D0                      !
      END DO                               !
                                           !
      MMP=0.0D0                            ! :::::::::::::::::::::::::::::::::::
      SMASS=0.0D0                          !
      countp=0D0                           !
      Pmaxp=0.0D0
      Pmax=0.0D0

C SWEEPS ALONG Y AXIS CONTROL VOLUMES:                                      
C...........................................................................!
          K=2                            !
          KM1=1                          !
          KP1=K+1                        !
                                         ! ON First Internal Vcv, k=2
          DYVK=DYV(K)                    ! Y size of Vcv
          SVNK=SVN(K)                    ! weight coeffs for V-eqn
          SVSK=SVS(K)                    !

          IF (Inerp.eq.1) Inerl=1        ! Inertia at entrance
                                         !
!	write(6,*) "befor calcv1",k ; call flush(6)
          CALL CALCV                     ! Solve V-eqn on k row -> V(j, k)
!	write(6,*) "after calcv1" ; call flush(6)
                                         ! 
          Inerl=Inerldum

      DO J=1, NXT                        ! Save pressure coeffs for extended
          DVV(J)=DV(J)                   !  land correction
      END DO                             ! ::::::::::::::::::::::::::::::::::

      K=1                                !

  700 KM1=K                              ! Solve:
      K=KM1+1                            !
      KP1=K+1                            !
      IF (K.EQ.NYI) KP1=NYI              !
      DYPK=DYP(K)                        ! Y size of U & P cvs

!	write(6,*) "befor calcu",k ; call flush(6)
      CALL CALCU                         ! Solve U-eqn on K row -> U(j, k)
!	write(6,*) "befor calcu" ; call flush(6)

      IF (K.EQ.NYI) goto 800             ! --> Update Vinlet

      KM1=K                              !
      K=KM1+1                            !
      KP1=K+1                            !
      IF (K.EQ.NYI) KP1=NYI               !
      DYVK=DYV(K)                        ! Y size at Vcv
      SVNK=SVN(K)                        ! Weight coeffs.
      SVSK=SVS(K)                        !

!	write(6,*) "befor calcv2" ; call flush(6)
      CALL CALCV                         ! Solve V-eqn on k row -> V(j, k)
!	write(6,*) "after calcv2" ; call flush(6)

      KP1=K                              !
      K=KM1                              ! Correct:
      KM1=K-1                            !

	
!	write(6,*) "befor calcp" ; call flush(6)
      CALL CALCP(Pmax,Countp)            ! P(j,k), U(j,k), V(j,k), V(j,k+1)
!	write(6,*) "after calcp" ; call flush(6)
	call flush(6)

      GOTO 700                           !
C                                        ! ::::::::::::::::::::::::::::::::::
C :::::::::::::::::::::::::::::::::::::: !
                                         !
                                         !
 800  SVNK=SVN(2)                        !
      SVSK=SVS(2)                        ! Calculate Inlet Velocity &  
                                         ! Modify entrance Pressure
!	write(6,*) "befor edgev" ; call flush(6)
      CALL EDGEV(Pmax,DYP(1),SVNK,SVSK,INERP,Countp)
!	write(6,*) "after edgev" ; call flush(6)

c   !....................................!................................!
c                                        ! Modify P(NYI) for end seal cond.
      IF (Cseal.ne.0.0D0)  THEN		   ! TS 07/2000 Possible Cseal<0
        CALL ENDS(DYP(Nyi),Cdseal,alfp, Countp)
      END IF
C----------------------------------------!
                                         !
!	write(6,*) "befor flows" ; call flush(6)
      CALL FLOWS                         ! Calculates flows on recesses
!	write(6,*) "after flows" ; call flush(6)

C    !...........................................!
      IF (SVERB.EQ.1) THEN                       !
          WRITE (6, 1200) ITER, COuntp, Pmaxp, Pmax, SMASS    
      ELSE IF (SVERB.EQ.3) THEN                  !
          DUMYI=0                                !
!	write(6,*) "befor printuvp" ; call flush(6)
          CALL PRINTUVP(DUMYI, 1)                !
!	write(6,*) "after printuvp" ; call flush(6)
      END IF                                     !
      IF (DEVICE) THEN                           !
          IF (DVERB.EQ.1) THEN                   !
              WRITE (6, 1200) ITER, Countp, Pmaxp, Pmax,SMASS  
          ELSE IF (DVERB.EQ.3) THEN              !
!	write(6,*) "befor printuvp" ; call flush(6)
              CALL PRINTUVP(DEVICE, 2)           !
!	write(6,*) "after printuvp" ; call flush(6)
          END IF                                 !
      END IF                                     !
C    !...........................................!

      IF (QIN.EQ.0.0D0) THEN
         SMASS=SMASS
      ELSE
         SMASS=SMASS/DABS(QIN)              ! Absolute mass flow error
      END IF

      if (iconv.ge.3) then               !
        if((Smass-Sold).gt.(0.0)) goto 1570    !-> error grows
      else                               !
        if(iter.eq.1) goto 100           ! 2 its minimum
        iconv=iconv+1                    !..........
      end if                             !
      Sold=Smass                         !..........
                                         !
	call flush(6)
C    !...................................!......................
                                         !
  900 continue                           ! Convergence in lands
C    !...................................!......................
      IF (countp.ge.Maxcountp) goto 1100 !
c    !...................................!.......................
                                         !
      IF (ITER.GE.ITMAX) goto 1000       ! not converged
C    !...................................!...................... 


      GOTO 100 ! Start                   !
                                         ! NEW iteration

C    !...................................!


C    !...................................!

 1570 alfp=alfp**1.25                    ! Reduce Relax-factor
      alfu=alfu**1.25                    !
      betu=1.0-Alfu                      !
      iconv=0                            !
      write(6,*)  'Changed relax Params: Ap & Au:', Alfp, Alfu,
     +             ' on iter:', iter
      IF ((Alfp.lt.(0.25)).or.(Alfu.lt.(0.25))) THEN
        write(6,*) 'Relax parameters too low: PRG Aborts in it:',iter
        goto 1000
      END IF                             !.............................!
      iter=iter-1
      goto 100                           ! -> START OVER

c........................................!.............................!

 1000 WRITE (6, *) 'WARNING:  MAXIMUM NUMBER OF ITERATIONS ON LANDS EXCE
     +EDED., ITMAX:', ITMAX 
      Pmaxp=Sflow
      alfp=relaxp
      alfu=relaxu
      betu=1.0-Alfu                      
      CALL BEEPER
      RETURN                             

                                          
 1100 Pmaxp=Sflow
      alfp=relaxp
      alfu=relaxu
      return  

 1200 FORMAT (' ', 3X, 'It:', I3, '  Countp:', I5,
     + 'Pmaxp:', E12.5E2, ' PMAX:', E12.5E2, 
     +                 '   SUM(Msourc):', E12.5E2)

 1300 FORMAT (' ', 'Convergence achieved in ', I3, ' iterations.')

      END 


C *****************************************************************************
C **                                                                         **
C **  Subroutine Setlim                                                      **
C **                                                                         **
C **  SETLIM:  Sets lower and upper limits for control volumes.              **
C **                                                                         **
C *****************************************************************************

      SUBROUTINE SETLIM

      IMPLICIT NONE

C ----------------------------------------------------------------------------
C --  COMMON variable declarations                                         --
C ----------------------------------------------------------------------------

      COMMON /UVEC/ JUMIN, JUMAX, JUSTART, JUSTOP
      COMMON /VVEC/ JVMIN, JVMAX, JVSTART, JVSTOP
      COMMON /PRES/ JPMIN, JPMAX, JPSTART, JPSTOP
      common /nodes/ nxt, nyi

      INTEGER JUMIN, JUMAX, JUSTART, JUSTOP
      INTEGER JVMIN, JVMAX, JVSTART, JVSTOP                                
      INTEGER JPMIN, JPMAX, JPSTART, JPSTOP
      integer nxt, nyi


C ----------------------------------------------------------------------------
C --  SETLIM code                                                           --
C ----------------------------------------------------------------------------

C                                 !.................................
                                  ! Inter=0, on 2PI extended lands
C                                 !.................................
          JUMIN=1                 ! min: leftmost point 
          JUMAX=NXT-1             ! max: rightmost point
          JUSTART=NXT-1           ! start: loc. left BC
          JUSTOP=NXT              ! stop: loc. right BC
          JPMIN=JUMIN             !
          JPMAX=JUMAX             !
          JPSTART=JUSTART         !
          JPSTOP=JUSTOP           !
          JVMIN=1                 !
          JVSTART=NXT-1           !
          JVMAX=JUMAX             !
          JVSTOP=JUSTOP           !
C                                 !................................
      
      END 


     
C *****************************************************************************
C **                                                                         **
C **  Subroutine Calcu                                                       **
C **                                                                         **
C **  CALCU:  Finds U velocity from marching equation:                       **
C **    Apu Up = Aep Ue + Awu Uw + Asu Us + Dpu (Pp-Pe) + Spu + (1-a)Apu Up  **
C **                                                                         **
C *****************************************************************************

      SUBROUTINE CALCU

      IMPLICIT NONE

      INCLUDE 'params.f'

C ----------------------------------------------------------------------------
C --  COMMON variable declarations                                         --
C ----------------------------------------------------------------------------
                                                                       
      COMMON /DXVEC/ DXP(MAXNXT), DXU(MAXNXT), SUW(MAXNXT),SUE(MAXNXT)
      COMMON /HFILM/ HP(MAXNXT,MAXNYI), HU(MAXNXT,MAXNYI),
     +               HV(MAXNXT,MAXNYI)
      COMMON /UVARRAY/ U(MAXNXT,MAXNYI), V(MAXNXT,MAXNYI)
      COMMON /PARRAY/ P(MAXNXT,0:MAXNYI)
      COMMON /RHOEMU/ RHP(MAXNXT,0:MAXNYI),EMP(MAXNXT,0:MAXNYI)
      COMMON /DPUV/ BU(MAXNXTP2), BV(MAXNXTP2), BVV(MAXNXTP2)
      COMMON /TDMA0/ A(MAXNXTP2), B(MAXNXTP2), C(MAXNXTP2), D(MAXNXTP2)
      COMMON /FACTORS/ REP, REY, MSPEED, SPEED,
     +                 klosy, SWIRL, PCAV, ALFU, BETU, ALFP
      COMMON /KVALA/ DYVK, DYPK, SVNK, SVSK

      COMMON /UVEC/ JUMIN, JUMAX, JUSTART, JUSTOP
      COMMON /FLAGS/ INERL, INERP, ITURB, ICAV
      COMMON /NODES/ Nxt, Nyi
      COMMON /KVALB/ K, KM1, KP1
c........................................................................

      DOUBLE PRECISION DXP, DXU,SUW, SUE, HP, HU, HV,
     +                 U, V, P, RHP, EMP,
     +                 BU, BV, BVV,A, B, C, D,
     +                 REP, REY, MSPEED, SPEED,
     +                 klosy, SWIRL, PCAV, ALFU, BETU, ALFP,
     +                 DYVK, DYPK, SVNK, SVSK

      INTEGER JUMIN, JUMAX, JUSTART, JUSTOP,
     +        INERL, INERP, ITURB, ICAV,
     +        Nxt, Nyi, K, KM1, KP1


C ----------------------------------------------------------------------------
C --  Local variable declarations                                          --
C ----------------------------------------------------------------------------
      DOUBLE PRECISION DYPRE, SANB, UW, UE, FWU, AEU,
     +                 VS, VN, VP, KXP, KXS, SPU, SPS, DPU, FEU,
     +                 FSU, FNU, ANU, AWU, ASU,
     +                 FNKX, ZERO,  US, Hs, Hn, Pp,
     +                 RHOP,EMUP,RHON,EMUN,RHOS,EMUS,RHOW,EMUW,
     +                 RHOE,EMUE,RHSE,EMSE,RHNE,EMNE,SUEJ,SUWJ
      INTEGER J, JJ, JJJ, JJMAX, JM1, JP1

C ----------------------------------------------------------------------------
C--  CALCU code                                                            --
C ----------------------------------------------------------------------------
C                                                !
      JJJ=2-JUMIN                                !
      JJMAX=JUMAX+JJJ

      DYPRE=DYPK*REY                             ! Dypk: Size of Ucv

      J=JUMIN 
      ZERO=0.0D0
      

                                                 ! BCs On extended lands:
c     !....................!                     !....
          UW=U(JUSTART, K)                       ! Left U velocity=BC
          UE=U(JUSTOP, K)                        ! Right U velocity=BC
          Rhow=Rhp(jumin,k)                      !
          FWU=DYPRE*RHOW*HP(JUMIN,K)*(UW+U(JUMIN,K))/2.0D0    ! west flow on Ucv
C     !.....................!                    ! 


     
      do j=jumin, jumax                          !........................
c    !......................!                    ! ON Ucv
 800  JP1=J+1                                    ! Sweep from left -> right
      JJ=J+JJJ                                   ! J+2-JUMIN
      SUEJ=SUE(J)                                !
      SUWJ=SUW(J)                                !
      
      VS=V(J,K)*SUEJ+V(JP1,K)*SUWJ               ! South vel=Vp for upwind
      VN=V(J,KP1)*SUEJ+V(JP1,KP1)*SUWJ           ! North vel=Vp for upwind
      VP=(VS+VN)/2.0D0                           !

      Rhop=Suej*Rhp(j,k)+Suwj*Rhp(jp1,k)
      Emup=Suej*Emp(j,k)+Suwj*Emp(jp1,k)

      KXP=FNKX(U(J,K),VP,HU(J,K),SPEED,REP,RHOP,EMUP,KXS)

      SPU=KXP*DXU(J)*DYPK/HU(J,K)                ! Source term, friction
      SPS=(KXS/KXP)*SPU                          !
      B(JJ)=SPU/ALFU                             ! Spu/Alfu
      C(JJ)=ZERO                                 !
      A(JJ)=ZERO                                 !
      DPU=HU(J,K)*DYPK                           !
      D(JJ)=DPU*(P(J,K)-P(JP1,K))+SPS*MSPEED     ! = Dpv (Ps-Pa)+U/2 Sps
      BU(J)=B(JJ) 


C     !....................!                     ! Account for Inertia
      IF (INERL.EQ.1) THEN                       ! Effects on lands: 
C     !....................!

          AWU=DMAX1(+FWU,ZERO)                   !At west Ucv face
       
          Rhoe=Rhp(jp1,k)                        !At east Ucv face
          FEU=DYPRE*RHOE*HP(JP1,K)*(U(J,K)+U(JP1,K))/2.0D0
          AEU=DMAX1(-FEU,ZERO)                   !


          HS=HV(J,K)*SUEJ+HV(JP1,K)*SUWJ         !At south Ucv face
          RHOS=(Rhp(j,km1)+Rhp(j,k))/2.0D0
          RHSE=(Rhp(jp1,km1)+Rhp(jp1,k))/2.0D0
          RHOS=RHOS*SUEJ+RHSE*SUWJ               !
          FSU=REY*RHOS*DXU(J)*HS*VS              !
          ASU=DMAX1(+FSU, ZERO)                  !

          HN=HV(J,KP1)*SUEJ+HV(JP1,KP1)*SUWJ     !At north Ucv face
          RHON=(Rhp(j,k)+Rhp(j,kp1))/2.0D0
          RHNE=(Rhp(jp1,k)+Rhp(jp1,kp1))/2.0D0
          RHON=RHON*SUEJ+RHNE*SUWJ               !
          FNU=REY*RHON*DXU(J)*HN*VN              !
          ANU=DMAX1(-FNU, ZERO)                  !            

          SANB=AEU+AWU+ASU+ANU                   ! SUM(Anb)
          C(JJ)=-AWU                             ! Coef. of Uw
          A(JJ)=-AEU                             ! Coef. of Ue
          B(JJ)=B(JJ)+SANB/ALFU                  ! Coef. of Up = Apu/Alfu

          US=0.5D0*(U(J, KM1)+U(J, K))           !
          D(JJ)=D(JJ)+ASU*US+ANU*U(J,KP1)        ! RHS coef.
          FWU=FEU                                !
          BU(J)=B(JJ)-SANB                       ! for SIMPLEC corection

      END IF                                     ! Inerl=1
C    !......................!                    !                           !

      D(JJ)=D(JJ)+BETU*B(JJ)*U(J, K)             ! RHS+Underrelax, Betu=1-Alfu

c    !......................!                    !
      END DO                                     ! J=jumin , jumax      
C  ..............................................!.......................
                                                 !
C                                                ! ::::::::::::::::::::::::::
C .............................................. ! Solution of U-equations by
  500 CALL TDMA(UW, UE, JJMAX)                   ! TDMA algorithm
C .............................................. ! ..........................
C                                                !
  600 DO J=JUMIN, JUMAX                          ! B(jj)=Usol(j)
          JJ=J+JJJ                               !
          U(J, K)=B(JJ)                          ! Order solution & set
          BU(J)=HU(J,K)*DYPK/BU(J)               ! coeffics. for Peqn.
      END DO                                     ! ..........................
C                                                !
C .............................................. ! .......................... 
          U(JUSTOP, K)=U(JUMIN, K)               ! periodicity condition
c
c................................................!...........................
      END 


C *****************************************************************************
C **                                                                         **
C **  Subroutine Calcv                                                       **
C **                                                                         **
C **  CALCV:  Finds V velocity from marching equation:                       **
C **       Apu Vp = Aep Ve + Awu Vw + Asu Vs + Dpu (Ps-Pp) + (1-a) Apu Vp*   **
C **                                                                         **
C *****************************************************************************

      SUBROUTINE CALCV

      IMPLICIT NONE

      INCLUDE 'params.f'


C ----------------------------------------------------------------------------
C --  COMMON variable declarations                                         --
C ----------------------------------------------------------------------------

      COMMON /DXVEC/ DXP(MAXNXT),DXU(MAXNXT),SUW(MAXNXT),SUE(MAXNXT) 
      COMMON /HFILM/ HP(MAXNXT,MAXNYI), HU(MAXNXT,MAXNYI),
     +               HV(MAXNXT,MAXNYI) 
      COMMON /UVARRAY/ U(MAXNXT,MAXNYI), V(MAXNXT,MAXNYI)
      COMMON /PARRAY/ P(MAXNXT,0:MAXNYI)
      COMMON /RHOEMU/ RHP(MAXNXT,0:MAXNYI), EMP(MAXNXT,0:MAXNYI)
      COMMON /DPUV/ BU(MAXNXTP2), BV(MAXNXTP2), BVV(MAXNXTP2)
      COMMON /TDMA0/ A(MAXNXTP2), B(MAXNXTP2), C(MAXNXTP2), D(MAXNXTP2)

      COMMON /FACTORS/ REP, REY, MSPEED, SPEED,
     +                 klosy, SWIRL, PCAV, ALFV, BETV, ALFP
      COMMON /KVALA/ DYVK, DYPK, SVNK, SVSK

      COMMON /VVEC/ JVMIN, JVMAX, JVSTART, JVSTOP
      COMMON /FLAGS/ INERL, INERP, ITURB, ICAV
      COMMON /NODES/ Nxt, Nyi
      COMMON /KVALB/ K, KM1, KP1

      DOUBLE PRECISION DXP, DXU, SUW, SUE, HP, HU, HV,
     +                 U, V, P, RHP, EMP,
     +                 BU, BV, BVV,A, B, C, D,
     +                 REP, REY, MSPEED, SPEED,
     +                 klosy, SWIRL, PCAV, ALFV, BETV, ALFP,
     +                 DYVK, DYPK, SVNK, SVSK

      INTEGER JVMIN, JVMAX, JVSTART, JVSTOP,
     +        INERL, INERP, ITURB, ICAV,
     +        Nxt, Nyi, K, KM1, KP1

C ----------------------------------------------------------------------------
C --  Local variable declarations                                          --
C ----------------------------------------------------------------------------
      DOUBLE PRECISION DYVRE, SANB, VW, VE, UW, FWV, UE, UP,
     +                 KYP, SPV, DPV, FEV, FSV, AEV, AWV, ASV,
     +                 ANV, c1, c2, FEW, VS, FSVO, ASVO,
     +                 CORRS, FORCEP, ZERO, FNV, FACTOR,
     +                 RHOP,EMUP,RHOS,RHON,RHOW,
     +                 RHSW,RHSE,RHOE,FNKY, HW,HE, AY,VN
                      
      INTEGER J, JJ, JJJ, JJMAX, JP1
	integer imrf
	imrf=0
C ----------------------------------------------------------------------------
C --  CALCV code                                                            --
C ----------------------------------------------------------------------------
!	write(6,*) "in calcv 1",imrf ; imrf=imrf+1; call flush(6)
      JJJ=2-JVMIN                                 !
      JJMAX=JVMAX+JJJ                             !
      DYVRE=DYVK*REY                              ! Dyvk = Size of Vcv
      ZERO=0.0D0                                  ! 

      IF (K.eq.2) THEN                            ! weight coefficients for
           c1=1.0D0                               ! Vsouth
           c2=0.0D0                               !
      else                                        !
           c1=0.5D0                               !
           c2=0.5D0                               !
      end if                                      !.......................


!	write(6,*) "in calcv 2",imrf ; imrf=imrf+1; call flush(6)
 100  VW=V(JVSTART, K)                            ! Left V velocity = BC
      VE=V(JVSTOP, K)                             ! right V velocity = BC
C                                                 !
      UW=U(JVSTART, K)*SVSK+U(JVSTART, KM1)*SVNK  ! West velocity on Vcv
      HW=HU(JVSTART,K)*SVSK+HU(JVSTART,KM1)*SVNK  !
      
      J=Jvstart
      RHOW=Rhp(j,k)*Sue(j)+Rhp(jvmin,k)*Suw(j)
      RHSW=Rhp(j,km1)*Sue(j)+Rhp(jvmin,km1)*Suw(j)
      RHOW=RHOW*SVSK+RHSW*SVNK                    !

      FWV=DYVRE*RHOW*HW*UW*INERL                  ! West flow on Vcv
C                                                 ! .........................
C ......................                          ! On standard V equation
!	write(6,*) "in calcv 3",imrf ; imrf=imrf+1; call flush(6)
      DO J=JVMIN, JVMAX                           ! .........................
C ......................                          !
          JJ=J+JJJ                                ! SWEEP from left -> right
          C(JJ)=ZERO                              ! on V control volume
          A(JJ)=ZERO                              !
          UE=U(J, K)*SVSK+U(J, KM1)*SVNK          ! East velocity
          UP=(UE+UW)/2.0D0                        ! U velocity at center Vcv

          RHOP=Rhp(j,k)*Svsk+Rhp(j,km1)*Svnk
          EMUP=Emp(j,k)*Svsk+Emp(j,km1)*Svnk

          KYP=FNKY(UP,V(J,K),HV(J,K),SPEED,REP,RHOP,EMUP)

          SPV=KYP*DXP(J)*DYVK/HV(J,K)             ! Source term, friction
          B(JJ)=SPV/ALFV                          ! Alfv=Underrelaxation par

          DPV=HV(J,K)*DXP(J)                      ! Area of pressure action
          D(JJ)=DPV*(P(J,KM1)-P(J,K))             ! = Dpv (Ps-Pp)

          BV(J)=B(JJ)
C                                                 !
C .............................                   ! Account for inertia
          IF (INERL.EQ.1) THEN                    ! Effects on lands:
C .............................                   !.................
                                                  !ON Vcv south face
              VS=C2*V(J,K)+C1*V(J,KM1)            !
              FSV=RHP(j,km1)*REY*DXP(J)*HP(J,KM1)*VS    
              ASV=DMAX1(+FSV,ZERO)                !
                                                  !......
              VN=0.5D0*(V(J,K)+V(J,KP1))          ! ON Vcv North face
              FNV=RHP(j,k)*REY*DXP(J)*HP(J,K)*VN      
              ANV=DMAX1(-FNV,ZERO)                !
                                                  !......
              HE=HU(J,K)*SVSK+HU(J,KM1)*SVNK      ! ON Vcv east face
              JP1=J+1                             !
              RHSE=Rhp(j,km1)*Sue(j)+Rhp(jp1,km1)*Suw(j)
              RHOE=Rhp(j,k)*Sue(j)+Rhp(jp1,k)*Suw(j)
              RHOE=RHOE*SVSK+RHSE*SVNK            !
              FEV=DYVRE*UE*HE*RHOE                !
              AEV=DMAX1(-FEV,ZERO)                !
                                                  !
              AWV=DMAX1(+FWV,ZERO)                !
    
              SANB=AEV+AWV+ASV+ANV                ! SUM(Anb)
              C(JJ)=-AWV                          ! coef. of Vw
              A(JJ)=-AEV                          ! coef. of Ve
              B(JJ)=B(JJ)+SANB/ALFV               ! coef. of Up = Apv/Alfv
              D(JJ)=D(JJ)+ASV*V(J,KM1)+ANV*V(J,KP1) ! RHS coef=Dpv(Ps-Pe)+Asv Vs

              FWV=FEV                             !
              BV(J)=B(JJ)-SANB                    ! For SIMPLEC correction
          END IF                                  ! Inerl=1
C .............................                   !

          D(JJ)=D(JJ)+BETV*B(JJ)*V(J, K)          ! RHS+Underrelax, Betv=1-Alfv

          UW=UE                                   !
C .............................                   !
      END DO                                      ! :::::::::::::::::::::::::
!	write(6,*) "in calcv 4",imrf ; imrf=imrf+1; call flush(6)
C .............................                   !

  200 IF ((INERP+INERL).EQ.0) GOTO 300            ! -> SOLVE equations
      IF (K.GT.2) GOTO 300                        !

 	
      DO J=JVMIN, JVMAX                           ! Modify eqns. at entrance
c   !..........................!                  !..........................

              JJ=J+JJJ                            ! Correct Flow & Coefs.  
              Vs=V(j,1)                           !

              Rhos=RHP(j,1)                       !    
              FSV=Rey*RHOS*Dxp(J)*HP(J,1)*Vs      ! south flow
              FSVo=FSV*Inerl                      !
              ASV=DMAX1(FSV,ZERO)                 !
              ASVo=DMAX1(Fsvo,ZERO)               !
              CORRS=(ASV-ASVO)/ALFV               ! Correction on south coefs.
              B(JJ)=B(JJ)+CORRS                   ! node
!	write(6,*) "in calcv 5",imrf ; imrf=imrf+1; call flush(6)
              BV(J)=BV(J)+CORRS*BETV              ! Simplec procedure
              FORCEP=DXP(J)*HV(J,K)*(P(J,1)-P(J,K))  ! = Dpv (Ps-Pp)
              D(JJ)=BETV*B(JJ)*V(J,1)             ! RHS side = Under-relax V

         IF ((INERP.eq.1).AND.(VS.GT.ZERO)) THEN  !............. Correct:
                  AY=Klosy*Rhos                   ! 
                  FACTOR=AY*DXP(J)*HV(J,K)*VS     !
                  ASV=ASV-FACTOR                  !  
                  BV(J)=BV(J)+FACTOR              !       
                  FORCEP=DXP(J)*HV(J,K)*(P(J,0)-P(J,K)) ! Dpv (Prec-Pp)
         END IF                                   ! Inerp=1
!	write(6,*) "in calcv 6",imrf ; imrf=imrf+1; call flush(6)

              VN=(V(J,K)+V(J,KP1))/2.0D0          !
              FNV=REY*RHP(j,k)*DXP(J)*HP(J,K)*VN*INERL
              ANV=DMAX1(-FNV,ZERO)                !
              D(JJ)=D(JJ)+ASV*VS+FORCEP+ANV*V(J,KP1)

          END DO                                  ! 
C      !.......................!                  !------
      
!	write(6,*) "in calcv 7",imrf ; imrf=imrf+1; call flush(6)


C                                                 !
C ............................................... ! Solution of V-equations by
  300 CALL TDMA(VW, VE, JJMAX)                    ! TDMA algorithm
C ............................................... ! .........................
  400 DO J=JVMIN, JVMAX                           ! Order V solution
          JJ=J+JJJ                                !
          V(J, K)=B(JJ)                           !
          BV(J)=HV(J,K)*DXP(J)/BV(J)              !
      END DO                                      !
!	write(6,*) "in calcv 8",imrf ; imrf=imrf+1; call flush(6)

      V(JVSTOP,K)=V(JVMIN,K)                      ! periodicity condition

!	write(6,*) "leaving calcv 9",k ; call flush(6)
      END 

C *****************************************************************************
C **                                                                         **
C **  Subroutine Calcp                                                       **
C **                                                                         **
C **  CALCP:  Solves Pres Correc Eqn at inter-recess lands, Inter=1:         **
C **          Ap Pp = Ae Pe' + Aw Pw' + As Ps' + Mp                          **
C **                                                                         **
C *****************************************************************************

      SUBROUTINE CALCP(Pmax,Countp)

      IMPLICIT NONE

      INCLUDE 'params.f'

C ----------------------------------------------------------------------------
C --  COMMON variable declarations                                         --
C ----------------------------------------------------------------------------

      COMMON /DXVEC/ DXP(MAXNXT), DXU(MAXNXT),SUW(MAXNXT),SUE(MAXNXT) 
      COMMON /DYVEC/ DYP(MAXNYI), DYV(MAXNYI),SVN(MAXNYI),SVS(MAXNYI)
      COMMON /HFILM/ HP(MAXNXT,MAXNYI), HU(MAXNXT,MAXNYI),
     +               HV(MAXNXT,MAXNYI) 
      COMMON /UVARRAY/ U(MAXNXT,MAXNYI), V(MAXNXT,MAXNYI)
      COMMON /PARRAY/ P(MAXNXT,0:MAXNYI)
      COMMON /RHOEMU/ RHP(MAXNXT,0:MAXNYI), EMP(MAXNXT,0:MAXNYI)
      COMMON /DPUV/ DU(MAXNXTP2), DV(MAXNXTP2), DVV(MAXNXTP2)
      COMMON /TDMA0/ A(MAXNXTP2), B(MAXNXTP2), C(MAXNXTP2), D(MAXNXTP2)

      COMMON /FACTORS/ REP, REY, MSPEED, SPEED,
     +                 klosy, SWIRL, PCAV, ALFU, BETU, ALFP
      COMMON /KVALA/ DYVK, DYPK, SVNK, SVSK
      COMMON /SOURCEA/ SMASS,PEPS, MMP, Pmaxp 

      COMMON /PRES/ JPMIN, JPMAX, JPSTART, JPSTOP
      COMMON /FLAGS/ INERL, INERP, ITURB, ICAV
      COMMON /KVALB/ K, KM1, KP1

      DOUBLE PRECISION DXP, DXU, HP, HU, HV, SUW, SUE,
     +                 DYP,DYV,SVN,SVS,
     +                 U, V, P, RHP, EMP,
     +                 DU, DV, DVV,A, B, C, D
      DOUBLE PRECISION REP, REY, MSPEED, SPEED,
     +                 klosy, SWIRL, PCAV, ALFU, BETU, ALFP,
     +                 DYVK, DYPK, SVNK, SVSK,
     +                 SMASS, PEPS, MMP, Pmaxp 

      INTEGER JPMIN, JPMAX, JPSTART, JPSTOP,
     +       INERL, INERP, ITURB, ICAV, K, KM1, KP1
       

C ----------------------------------------------------------------------------
C --  Local variable declarations                                          --
C ----------------------------------------------------------------------------

      INTEGER Countp, J, JJ, JJJ, JJMAX
      DOUBLE PRECISION DPU, FWP, AWP, DPV, FEP, FSP, FNP, MP,
     +                 AEP, ASP, ANP, PW, PE, PPJK, FSNP,
     +                 RHOW,RHOE,RHOS,RHON,
     +                 Pmax, SVNKM1,SVSKM1, DMP, Ppmag,Pabs

C ----------------------------------------------------------------------------
C --  CALCP code                                                            --
C ----------------------------------------------------------------------------
      JJJ=2-JPMIN                              !
      JJMAX=JPMAX+JJJ                          !
      JJ=JPMIN+JJJ                             !
      J=JPSTART                                !

      RHOW=Rhp(j,k)*Sue(j)+Rhp(jpmin,k)*Suw(j)
      DPU=HU(J,K)*DYPK*RHOW                    !
      FWP=DPU*U(J, K)                          ! Flow on west side of Pcv
      AWP=0.0D0                                ! Awp coeff.


      KM1=K-1                                  !
      SVNKM1=SVN(K)                            !
      SVSKM1=SVS(K)                            ! 

C ......................                       !
      DO J=JPMIN, JPMAX                        !----------------------------
C ......................                       !
          JJ=J+JJJ                             ! Sweep from left-right on
                                               ! P cvs

          Rhoe=Rhp(j,k)*Sue(j)+Rhp(j+1,k)*Suw(j)
          DPU=HU(J,K)*DYPK*RHOE                ! east U flow on Pcv
          FEP=DPU*U(J,K)                       !
          AEP=DPU*DU(J)                        !

          Rhon=Rhp(j,kp1)*Svsk+Rhp(j,k)*Svnk   !
          DPV=HV(J,KP1)*DXP(J)*RHON            !
          FNP=DPV*V(J,KP1)                     ! North V flow on Pcv
          ANP=DPV*DV(J)                        !

          Rhos=Rhp(j,k)*Svskm1+Rhp(j,km1)*Svnkm1
          DPV=HV(J,K)*DXP(J)*RHOS              !
          FSP=DPV*V(J,K)                       ! South V flow on Pcv
          ASP=DPV*DVV(J)                       !

          MP=FWP-FEP+FSP-FNP                   ! Mass source on Pcv

          C(JJ)=-AWP                           ! coeffs. for TDMA soln.
          A(JJ)=-AEP                           !
          B(JJ)=AEP+AWP+ASP+ANP                !
          D(JJ)=MP                             ! P'(j, k-1)=0

          FWP=FEP                              !
          AWP=AEP                              !
          DMP=DABS(MP)                         ! 
          SMASS=SMASS+DMP                      ! Global mass source
          IF (DMP.GT.MMP) MMP=DMP              ! MAX(Local Mass source)
C
C......................                        !

      END DO                                   !
C ......................                       ! ::::::::::::::::::::::::::::
C                                              !

          B(JJMAX)=B(JJMAX)+A(JJMAX)           ! On extended lands, Ue spec.
          A(JJMAX)=0.0D0                       ! & Aep=0
                                               !
C                                              ! ............................
  100 PW=0.0D0 !=Pcor(Jpstart)                 ! Spec. pressure at left B
      PE=0.0D0 !=Pcor(Jpstop)                  ! Spec. pressure at right B
C                                              ! ----------------------------
  200 CALL TDMA(PW, PE, JJMAX)                 ! Solve by TDMA algorithm
C                                              ! ----------------------------
          B(JJMAX+1)=B(2)                      ! P'(Nxt)=P'(1)
                                               ! where Pcor(j)=B(jj)
C ............................................ !
C                                              ! ::::::::::::::::::::::::::::
  300 DO J=JPMIN, JPMAX                        ! CORRECT:
C    !...................!                     !
          JJ=J+JJJ                             !
          PPJK=B(JJ)                           ! P'(j, k)
          Ppmag=DABS(Ppjk)                     !

          IF (Ppmag.le.Peps) Countp=Countp+1

          IF (Ppmag.gt.Pmaxp) Pmaxp=Ppmag

          P(J, K)=P(J, K)+ALFP*PPJK            ! Corrected P(j, k),

          Pabs=DABS(P(j,k))

          IF (Pabs.gt.Pmax) Pmax=Pabs          !

          U(J, K)=U(J, K)+(PPJK-B(JJ+1))*DU(J) ! Corrected U(j, k) velocity
          V(J, K)=(V(J, K)-PPJK*DVV(J))        !     "     south V velocity
          V(J, KP1)=(V(J, KP1)+PPJK*DV(J))     !     "     north V velocity
          DVV(J)=DV(J)                         ! for use on next k
 

      END DO                                   !
C    !...................!                     !
C                                              !
C ............................................ ! ----------------------------
C                                              !
          P(JPSTOP, K)=P(JPMIN, K)             ! On extended lands
          V(JPSTOP, K)=V(JPMIN, K)             ! jpstop=nxt, jpmin=1
          U(JPSTOP, K)=U(JPMIN, K)             ! Inforce periodicity
          V(JPSTOP, KP1)=V(JPMIN, KP1)         !
c .............................................!.............................

  400 DO J=Jpmin, Jpmax                        ! Calculate Props
        Call Locprops(Rhp(j,k),Emp(j,k),P(j,k))
      END DO
        Rhp(jpstop,k)=Rhp(jpmin,k)
        Emp(jpstop,k)=Emp(jpmin,k)

      END 



C *****************************************************************************
C **                                                                         **
C **  Subroutine Edgev                                                       **
C **                                                                         **
C **  EDGEV:  Updates V velocity and entrance Pres  due to inertia effects.  **
C **                                                                         **
C *****************************************************************************

      SUBROUTINE EDGEV(Pmax,DYPK, SVNK, SVSK, INERP, Countp)

      IMPLICIT NONE

      INCLUDE 'params.f'


C ----------------------------------------------------------------------------
C --  COMMON variable declarations                                         --
C ----------------------------------------------------------------------------

      COMMON /DXVEC/ DXP(MAXNXT), DXU(MAXNXT),SUW(MAXNXT),SUE(MAXNXT) 
      COMMON /HFILM/ HP(MAXNXT,MAXNYI), HU(MAXNXT,MAXNYI),
     +               HV(MAXNXT,MAXNYI)
      COMMON /UVARRAY/ U(MAXNXT,MAXNYI), V(MAXNXT,MAXNYI)
      COMMON /PARRAY/ P(MAXNXT,0:MAXNYI)
      COMMON /RHOEMU/ RHP(MAXNXT,0:MAXNYI), EMP(MAXNXT,0:MAXNYI)
      COMMON /FACTORS/ REP, REY, MSPEED, SPEED,
     +                 klosy, SWIRL, PCAV, ALFV, BETV, ALFP
      COMMON /SOURCEA/ SMASS, PEPS, MMp, PMAXP
      COMMON /NODES/ Nxt, Nyi

      DOUBLE PRECISION DXP, DXU, SUW, SUE, HP, HU, HV,
     +                 U, V, P, RHP, EMP,
     +                 REP, REY, MSPEED, SPEED,
     +                 klosy, SWIRL, PCAV, ALFV, BETV, ALFP,
     +                 SMASS, PEPS, MMp, PMAXP
      INTEGER Nxt, Nyi
C ----------------------------------------------------------------------------
C --  Local variable declarations                                          --
C ----------------------------------------------------------------------------

      DOUBLE PRECISION DYPK, SVSK, SVNK,
     +                 VS,Fsv, Fnv, Few, PEDGE,
     +                 Pmax, Pcorr, AY, RHOS,EMUS,RHON,EMUN
      INTEGER Countp, INERP, j, jm1

C ----------------------------------------------------------------------------
C --  EDGEV code                                                            --
C ----------------------------------------------------------------------------

C    !.........................!               !..........................
      DO J=2, Nxt                              ! SWEEP from
C    !.........................!               ! left Vcv to right vcv
         Jm1=J-1

              Rhos=Rhp(j,1)
              Rhon=Rhp(j,2)*Svsk+Rhos*Svnk
              Fnv=Rhon*Dxp(J)*Hv(j,2)*V(j,2)   
              Few=Rhos*Dypk*Swirl*(Hu(j,1)-Hu(jm1,1))

              Fsv=Fnv+Few
              Vs=Fsv/Rhos/Dxp(j)/Hp(j,1)
              PEDGE=P(j,0)                       != Pexternal
              Pcorr=0.0D0

              IF ((VS.GT.(0.0D0)).AND.(Inerp.eq.1)) THEN
c            !........................................!

                 PEDGE=PEDGE-klosy*Rhos*VS*VS
                 PEDGE=DMAX1(PEDGE,PCAV)

                 IF (Pedge.gt.Pmax) Pmax=Pedge
  
                 Pcorr=DABS(Pedge-P(j,1))
                 IF (Pcorr.gt.Pmaxp) Pmaxp=Pcorr

              
              END IF
c           !..........................................!

         IF(Pcorr.le.peps) Countp=Countp+1

         V(J,1)=Vs

         P(J,1)=Pedge

         call locprops(rhp(j,1),emp(j,1), Pedge)

C    !........................!  
      END DO                                   ! ----------------------------
C    !........................!                !

C                                              !
  400 V(1,1) =V(NXT,1)                         ! 1st Vcv = Nxt(last) Vcv
      P(1,1)=P(NXT,1)                          ! set periodicity conditions
      V(1,2)=V(NXT,2)                          ! ----------------------------
      Rhp(1,1)=Rhp(Nxt,1)
      Emp(1,1)=Emp(Nxt,1)

      END


C *****************************************************************************
C **                                                                         **
C **  Subroutine Flows                                                       **
C **                                                                         **
C **  FLOWS:  Calculates inlet and outlet flow from annular seal             **
C **                                                                         **
C *****************************************************************************

      SUBROUTINE FLOWS

      IMPLICIT NONE

      INCLUDE 'params.f'

C ----------------------------------------------------------------------------
C --  COMMON variable declarations                                         --
C ----------------------------------------------------------------------------
      COMMON /DXVEC/ DXP(MAXNXT), DXU(MAXNXT),SUW(MAXNXT), SUE(MAXNXT)
      COMMON /HFILM/ HP(MAXNXT,MAXNYI), HU(MAXNXT,MAXNYI), 
     +               HV(MAXNXT,MAXNYI)
      COMMON /DYVEC/ DYP(MAXNYI), DYV(MAXNYI), SVN(MAXNYI), SVS(MAXNYI)
      COMMON /UVARRAY/ U(MAXNXT,MAXNYI), V(MAXNXT,MAXNYI)
      COMMON /PARRAY/ P(MAXNXT,0:MAXNYI)
      COMMON /RHOEMU/ RHP(MAXNXT,0:MAXNYI), EMP(MAXNXT,0:MAXNYI)
      COMMON /FLOW0/ QIN, QOUT, QFACTOR

      COMMON /NODES/ Nxt, Nyi

      DOUBLE PRECISION DXP, DXU,SUW, SUE, HP,HU, HV,
     +                 DYP, DYV, SVN, SVS,
     +                 U, V, P, RHP, EMP,
     +                 QIN, QOUT, QFACTOR
      INTEGER Nxt,Nyi

C ----------------------------------------------------------------------------
C --  Local variable declarations                                          --
C ----------------------------------------------------------------------------
      INTEGER  J, NYM1
      DOUBLE PRECISION SVNK

C ----------------------------------------------------------------------------
C -- FLOWS code                                                             --
C ----------------------------------------------------------------------------
  100 QIN=0.0D0
      QOUT=0.0D0 
      SVNK=SVN(Nyi)                      !
      Nym1=Nyi-1                         !

C    !.................!                 !
      DO J=1, Nxt-1                      ! Integrate in circ. direction
C    !.................!                 !

        Qin=Qin+dxp(j)*Rhp(j,1)*Hp(j,1)*V(j,1)!

        Qout=Qout+dxp(j)*(Rhp(j,nym1)*Svnk+Rhp(j,Nyi)*Svs(Nyi))
     +                    *Hv(j,Nyi)*V(j,Nyi)

C     !..................!
       END DO
C     !..................!               !

      END 

C *****************************************************************************
C **                                                                         **
C **  FORCE:  Integrates pressure to find forces.                            **
C **                                                                         **
C *****************************************************************************

      SUBROUTINE FORCE(DEVICE)

      IMPLICIT NONE

      INCLUDE 'params.f'


C ----------------------------------------------------------------------------
C --  COMMON variable declarations                                         --
C ----------------------------------------------------------------------------

      COMMON /XYVEC/ XP(MAXNXT), XU(MAXNXT),YP(MAXNYI), YV(MAXNYI) 
      COMMON /DXVEC/ DXP(MAXNXT), DXU(MAXNXT),SUW(MAXNXT),SUE(MAXNXT)
      COMMON /DYVEC/ DYP(MAXNYI), DYV(MAXNYI), SVN(MAXNYI), SVS(MAXNYI)
      COMMON /TRIGS/ COSXP(MAXNXT), SINXP(MAXNXT), 
     +               COSXU(MAXNXT), SINXU(MAXNXT)
      COMMON /PARRAY/ P(MAXNXT,0:MAXNYI)

      COMMON /PARAM1/ CLEAR, DIAM, LENGTH
      COMMON /ALIGNM/ AXO, AYO, ZO
      COMMON /FACTORS/ REP, REY, MSPEED, SPEED,
     +                 klosy, SWIRL, PCAV, ALFU, BETU, ALFP
      COMMON /FORCE0/ FFACTOR, FX, FY, TO, TOR
      COMMON /MOMENT0/ MFACTOR, MX,MY
      COMMON /NODES/ Nxt, Nyi

      DOUBLE PRECISION P, DXP, DXU, SUW, SUE, XP,XU,YP,YV,
     +                 DYP, DYV, SVN, SVS, AXO, AYO, ZO,
     +                 COSXP, SINXP, COSXU, SINXU
      DOUBLE PRECISION CLEAR, DIAM, LENGTH,
     +                 REP, REY, MSPEED, SPEED,
     +                 klosy, SWIRL, PCAV, ALFU, BETU, ALFP,
     +                 FFACTOR, FX, FY, TO, TOR, MFACTOR,MX,MY
      INTEGER Nxt, Nyi
 
C ----------------------------------------------------------------------------
C --  Local variable declarations                                          --
C ----------------------------------------------------------------------------

      INTEGER DEVICE, ICHECK, J, K, KM1
      DOUBLE PRECISION DYVK, PP, PCL, PSL, AREA, PCR, PSR, LD,
     +                 LOAD, ANGLE, YO, YY
      DOUBLE PRECISION YYP, PPM,PCLM,PSLM,PCRM,PSRM

C ----------------------------------------------------------------------------
C --  FORCE code                                                            --
C ----------------------------------------------------------------------------

      FX=0.0D0                          !
      FY=0.0D0                          !
      MX=0.0D0                          !
      MY=0.0D0                          !
      YO=ZO*2.0D0/DIAM                  ! PIVOT for moments
C    !...............                   !

      DO K=2, NYI                       ! On EXTENDED LANDS
C    !...............                   !------------------
          KM1=K-1                       !
          DYVK=DYV(K)                   ! Y size VCV
          YY=YP(K)-YO
          YYP=YP(KM1)-YO

          J=1                           ! Start of integration
          PP=DMAX1(PCAV,P(J,K))+DMAX1(PCAV,P(J,KM1))
          PPM=DMAX1(PCAV,P(J,K))*YY+
     +        DMAX1(PCAV,P(J,KM1))*YYP
          PCL=PP*COSXP(J)               ! left P*COS(Angle)
          PSL=PP*SINXP(J)               !  ''  P*SIN(Angle)
          PCLM=PPM*COSXP(J)
          PSLM=PPM*SINXP(J)

          DO J=2, NXT                   !
              AREA=DYVK*DXU(J-1)/4.0D0  ! Area of integration/4

              PP=DMAX1(PCAV,P(J,K))+DMAX1(PCAV,P(J,KM1))
              PCR=PP*COSXP(J)           ! right P*COS(Angle)
              PSR=PP*SINXP(J)           !  ''   P*SIN(Angle)
              FX=FX+(PCR+PCL)*AREA      ! X-force
              FY=FY+(PSR+PSL)*AREA      ! Y-force

              PPM=DMAX1(PCAV,P(J,K))*YY+
     +            DMAX1(PCAV,P(J,KM1))*YYP
              PCRM=PPM*COSXP(J)           ! right P*COS(Angle)
              PSRM=PPM*SINXP(J)           !  ''   P*SIN(Angle)
              MY=MY+(PCRM+PCLM)*AREA      ! Y-Moment
              MX=MX-(PSRM+PSLM)*AREA      ! X-Moment

              PCL=PCR                   !
              PSL=PSR                   !
              PCLM=PCRM
              PSLM=PSRM

          END DO                        ! j=2, Nxt
                                        !
      END DO                            ! k=2, Nyi
C    !....................!             !

C.......................................!.....................

  200 FX=FFACTOR*FX                     ! forces Fx&Fy in N
      FY=FFACTOR*FY                     !
      MX=MFACTOR*MX                     ! Moments Mx&My in N-m
      MY=MFACTOR*MY                     !
      ICHECK=0                          !
      CALL PRINTF(DEVICE,ICHECK)        ! Print forces

      END      

C
C 
C *****************************************************************************
C **                                                                         **
C **  Subroutine Torque                                                      **
C **                                                                         **
C **  TORQUE:  Calculates friction torque in film lands.                     **
C **                                                                         **
C *****************************************************************************

      SUBROUTINE TORQUE(DEVICE)

      IMPLICIT NONE

      INCLUDE 'params.f'

C ----------------------------------------------------------------------------
C --  COMMON variable declarations                                         --
C ----------------------------------------------------------------------------

      COMMON /DXVEC/ DXP(MAXNXT), DXU(MAXNXT), SUW(MAXNXT),SUE(MAXNXT)
      COMMON /HFILM/ HP(MAXNXT,MAXNYI), HU(MAXNXT,MAXNYI),
     +               HV(MAXNXT,MAXNYI)
      COMMON /DYVEC/ DYP(MAXNYI), DYV(MAXNYI), SVN(MAXNYI), SVS(MAXNYI)
      COMMON /PARRAY/ P(MAXNXT,0:MAXNYI)
      COMMON /UVARRAY/ U(MAXNXT,MAXNYI), V(MAXNXT,MAXNYI)
      COMMON /RHOEMU/ RHP(MAXNXT,0:MAXNYI), EMP(MAXNXT,0:MAXNYI)
      COMMON /FACTORS/ REP, REY, MSPEED, SPEED, 
     +                 klosy, SWIRL, PCAV, ALFU, BETU, ALFP
      COMMON /NODES/ Nxt, Nyi 
      COMMON /FORCE0/ FFACTOR, FX, FY, TO, TOR
      COMMON /VERB/ SVERB, DVERB, BEEP

      DOUBLE PRECISION DXP, DXU, SUW, SUE, HP, HU, HV,
     +                 DYP, DYV, SVN, SVS,
     +                 U, V, P, RHP, EMP,
     +                 REP, REY, MSPEED, SPEED,
     +                 klosy, SWIRL, PCAV, ALFU, BETU, ALFP,
     +                 FFACTOR, FX, FY, TO, TOR
      INTEGER Nxt, Nyi, SVERB, DVERB, BEEP

C ----------------------------------------------------------------------------
C --  Local variable declarations                                          --
C ----------------------------------------------------------------------------

      INTEGER DEVICE, J, K, KP1, JP1
      DOUBLE PRECISION VS, H, VN, VP, GAMA, PE, PP, DX
      DOUBLE PRECISION FNGAMA,TORQ, RHOP,EMUP
C ----------------------------------------------------------------------------
C --  TORQUE code                                                           --
C ----------------------------------------------------------------------------

      TOR=0.0D0
                                        !

      DO J=1, NXT-1                     ! ---> Torque on EXTENDED lands
C    !....................!             !
          JP1=J+1                       ! -----------------------------
          VS=V(J,1)*SUE(J)+V(JP1,1)*SUW(J)
          DX=DXU(J)                     !
                                        !
          DO K=1, NYI-1                 ! at interior UCV of extented land
              KP1=K+1                   !
              VN=V(J, KP1)*SUE(J)+V(JP1, KP1)*SUW(J)
              VP=(VN+VS)/2.0D0 
              H=HU(J,K)                 
              
              Rhop=Rhp(j,k)*Sue(j)+Rhp(jp1,k)*Suw(j)
              EMup=Emp(j,k)*SUe(j)+Emp(jp1,k)*Suw(j)

              GAMA=FNGAMA(U(J,K),VP,H,SPEED,REP,RHOP,EMUP) 
              PE=DMAX1(PCAV,P(JP1,K))   !
              PP=DMAX1(PCAV,P(J,K))     !
              TOR=TOR+(H*(PE-PP)/2.0D0/DX+GAMA)*DX*DYP(K)
              VS=VN                     !
          END DO                        !
       
          Rhop=Rhp(j,Nyi)
          Emup=Emp(j,Nyi)
          H=Hu(j,Nyi)             
          GAMA=FNGAMA(U(J,NYI),VS,H,SPEED,REP,RHOP,EMUP) ! At outlet (Exit) of bearing
          TOR=TOR+GAMA*DX*DYP(NYI)      ! PE=PP=0.0

C     !....................!            !
      END DO                            !
C    !..................................!......................

      TORQ=TO*TOR                       !
      WRITE (6, 100) TORQ               !
      IF (DEVICE) THEN                  !
          WRITE (1, 100) TORQ           !
      END IF                            !
C.......................................!..............................

  100 FORMAT (' |    Torque on Film Lands=', E12.5E2, ' N-m', 36X, '|', 
     +        /, ' +', 77('-'), '+')

  200 FORMAT (' ', '+', 77('-'), '+')

      END 

C
C TURBULENT SHEAR PARAMETERS MOVED TO turbcoefs.f on 6/8/92
c

C *****************************************************************************
C **                                                                         **
C **  Subroutine Tdma                                                        **
C **                                                                         **
C **  TDMA:  Solution of algebraic tridiagonal system of equations.          **
C **                                                                         **
C *****************************************************************************

      SUBROUTINE TDMA(ZLEFT, ZRIGHT, JMAX)


      IMPLICIT NONE

      INCLUDE 'params.f'

C ----------------------------------------------------------------------------
C --  COMMON variable declarations                                         --
C ----------------------------------------------------------------------------

      COMMON /TDMA0/ A(MAXNXTP2), B(MAXNXTP2), C(MAXNXTP2), D(MAXNXTP2)

      DOUBLE PRECISION A, B, C, D

C ----------------------------------------------------------------------------
C --  Local variable declarations                                          --
C ----------------------------------------------------------------------------

      DOUBLE PRECISION ZLEFT, ZRIGHT, DEN
      INTEGER JMAX, J, JM1

C ----------------------------------------------------------------------------
C --  TDMA code                                                             --
C ----------------------------------------------------------------------------

      D(1)=ZLEFT                        ! = F(1)
      A(1)=0.0                          ! = E(1)
      DO J=2, JMAX                      ! 1st sweep:
          JM1=J-1                       ! Find recursive
          DEN=1.0/(C(J)*A(JM1)+B(J))    !
          A(J)=-A(J)*DEN                ! = E(J)
          D(J)=(D(J)-C(J)*D(JM1))*DEN   ! = F(J)
      END DO                            !
      B(JMAX+1)=ZRIGHT                  ! Set BC at right
      DO J=JMAX, 1, -1                  ! 2nd sweep:
          B(J)=A(J)*B(J+1)+D(J)         ! calculate soln.
      END DO                            ! stored in vector B
      
      END 

C *****************************************************************************
C **                                                                         **
C **  Subroutine Ends                                                        **
C **                                                                         **
C **  ENDS   : Updates V velocity and discharge pressure due to end seal  .  **
C **                                                                         **
C *****************************************************************************

      SUBROUTINE ENDS(DYPK,Sealcoef,alfp, Countp)

      IMPLICIT NONE

      INCLUDE 'params.f'


C ----------------------------------------------------------------------------
C --  COMMON variable declarations                                         --
C ----------------------------------------------------------------------------
      COMMON /DXVEC/ DXP(MAXNXT), DXU(MAXNXT),SUW(MAXNXT),SUE(MAXNXT)
      COMMON /UVARRAY/ U(MAXNXT,MAXNYI), V(MAXNXT,MAXNYI)
      COMMON /PARRAY/  P(MAXNXT,0:MAXNYI)
      COMMON /RHOEMU/  RHP(MAXNXT,0:MAXNYI), EMP(MAXNXT,0:MAXNYI)
      COMMON /HFILM/ HP(MAXNXT,MAXNYI), HU(MAXNXT,MAXNYI),
     +               HV(MAXNXT,MAXNYI)
      COMMON /SOURCEA/ SMASS, PEPS, MMp, PMAXP
      COMMON /NODES/ NXT, NYI
c     .................................................................
      DOUBLE PRECISION DXP, DXU, SUW, SUE, HP, HU, HV,Pcorr
      DOUBLE PRECISION U,V,P,RHP,EMP,SMASS, PEPS, MMp, PMAXP
      INTEGER NXT, NYI, Countp

C ----------------------------------------------------------------------------
C --  Local variable declarations                                          --
C ----------------------------------------------------------------------------
      DOUBLE PRECISION zero,DYPK,  Psump, Sealcoef
      DOUBLE PRECISION Alfp, Rhoe,Rhow, Fe, Fw, Vend, Pend
      INTEGER J,KK

C ----------------------------------------------------------------------------
C --  ENDS CODE                                                            --
C ----------------------------------------------------------------------------
      KK=NYI
      zero=0.0D0
      Psump=zero               ! sump pressure

      j=NXT-1
      Rhow=Sue(j)*Rhp(j,kk)+Suw(j)*Rhp(1,kk)
      Fw=DYPK*Rhow*U(j,KK)*Hu(j,KK)
 
C    !.........................!               !..........................
      DO J=1, NXT-1                            ! Sweep from
C    !.........................!               ! left Vcv to right vcv
        Rhoe=Sue(j)*Rhp(j,kk)+Suw(j)*Rhp(j+1,kk)
        Fe=DYPK*Rhoe*U(j,KK)*Hu(j,KK)

c        Vend=(V(j,KK)*HV(j,KK)-(Fe-Fw)/Dxp(j)/Rhp(j,KK))/HP(j,KK)
         Vend=V(j,KK)

              IF (Vend.GT.zero) THEN
                 Pend=Psump+Sealcoef*Vend*Vend*Rhp(j,KK)
                 Pend=P(j,KK)*(1.0D0-Alfp)+Pend*Alfp
              ELSE
                 Pend=Psump
              END IF
                 Pcorr=DABS(Pend-P(j,KK))
                 P(j,KK)=Pend
                 IF (Pcorr.le.PEPS) Countp=Countp+1

         Fw=Fe

         CALL Locprops(Rhp(j,KK),Emp(j,KK),Pend)
C    !........................!
      END DO                                   ! ----------------------------
C    !........................!                !

      P(NXT,KK) =P(1,KK)                        ! set periodicity conditions
      Rhp(NXT,KK)=Rhp(1,KK)
      Emp(NXT,KK)=Emp(1,KK)

      END


C:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
C                                 
C END of CALCSOLNs.FOR for annular seals/ with end restriction
C
C
C:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
                                        
