      Program SIMCZO

      INTEGER          DIM,IDIV,JDIV,IDIM,JDIM,MXNOD,MXELE,NLOK,
     +                 IBMAX
      INTEGER          ICOORD, IELTOP

      PARAMETER       (DIM=2,IDIV=30,JDIV=50,IDIM=2*IDIV+1,
     +                 JDIM=2*JDIV+1,NLOK=8,IBMAX=12*IDIV+7)
      PARAMETER       (MXNOD=IDIM*JDIM,MXELE=IDIV*JDIV)
      PARAMETER       (ICOORD=DIM,IELTOP=10)

      DOUBLE PRECISION ZERO

      PARAMETER ( ZERO=0.0D0 )

      INTEGER  ELTOP(IELTOP,MXELE)
      INTEGER  STEER(10),NODEL 
      INTEGER  NS, J, NODE
      INTEGER  ITYP
      INTEGER  IP,INFILE

      DOUBLE PRECISION COORD(DIM,MXNOD),ELEMCO(DIM,NLOK)
      common /name1/ COORD, ELEMCO

      DOUBLE PRECISION A(2*MXNOD,IBMAX),F(2*MXNOD)
      common /name2/ A, F

      DOUBLE PRECISION ALOK(4*NLOK*NLOK),MLOK(4*NLOK*NLOK),
     +     FLOK(2*NLOK), GEOM(2*NLOK)
      common /name3/ ALOK, MLOK, FLOK, GEOM

      DOUBLE PRECISION U(2*MXNOD),T(MXNOD)

      DOUBLE PRECISION UPREV(2*MXNOD),TPREV(MXNOD),UELEM(2*NLOK),
     +		       TELEM(NLOK),FLOWEL(2*NLOK),FLX(NLOK),FLY(NLOK)
      common /name4/ U, T, UPREV, TPREV, UELEM, 
     +               TELEM, FLOWEL, FLX, FLY

      DOUBLE PRECISION VERR,TERR
      DOUBLE PRECISION EPS

      DOUBLE PRECISION RO,CP,DK(DIM,DIM),DLHEAT,XVEL,YVEL,
     +                 XLEN,YLEN,GRAVX,GRAVY,XP,YP,XO,YO,DT,
     +                 XVELO(NLOK),YVELO(NLOK),TGAMIN,TMELT,TGAMOU,
     +                 VISCOS,HEXPCO,REFTEM,RADCRY,RADCRU,ALFA,AREFT,
     +                 SURTEN,GASDEN,COANG,ASTEF,ALIQ,UEXT,DSIG

      DOUBLE PRECISION RELXM,RELXC,RELYM,RELYC,TIME

      DOUBLE PRECISION SINIT,SHCRU(IDIM),SHLIQ(IDIM),FREEB(IDIM),
     +                 FREEBO(IDIM)

      DOUBLE PRECISION ERRFB,TNOEPS

      INTEGER          KSW,NNODES,NELEMS,NELE,NNODX,NNODXC,NELEYM,
     +                 NELEYC,SIDNUM,ITCOUP,MAXCPL,
     +                 NELEXC,NELEXM,NELEMM,NNODM

      INTEGER          MAXIFB,ILIN1,ILIN2,INDMEL,IND1,IND2

      INTEGER          IBANDV,IDIAGV,LOV,IUPV,
     +                 IBANDT,IDIAGT,LOT,IUPT

      DOUBLE PRECISION ENORM,CONVT
      INTEGER NNN,IND,ITER
      DOUBLE PRECISION XXX,YYY

      OPEN (UNIT=7,FILE='param.dta',STATUS='OLD',ERR=990)
      READ(7,*) CP
      READ(7,*) DK(1,1),DK(1,2)
      READ(7,*) DK(2,1),DK(2,2)
      READ(7,*) 
      READ(7,*) SURTEN
      READ(7,*) GASDEN
      READ(7,*) COANG
      READ(7,*) 
      READ(7,*) HEXPCO
      READ(7,*) REFTEM
      READ(7,*) GRAVX
      READ(7,*) GRAVY
      READ(7,*) 
      READ(7,*) DLHEAT
      READ(7,*) 
      READ(7,*) XVEL
      READ(7,*) YVEL
      READ(7,*) 
      READ(7,*) TGAMIN
      READ(7,*) TMELT
      READ(7,*) TGAMOU
      READ(7,*) ALFA
      READ(7,*) AREFT
      READ(7,*) 
      READ(7,*) RADCRY
      READ(7,*) RADCRU
      READ(7,*) YLEN
      READ(7,*) SINIT
      READ(7,*) 
      READ(7,*) NELEXC,NELEXM
      READ(7,*) RELXC,RELXM
      READ(7,*) 
      READ(7,*) NELEYM,NELEYC
      READ(7,*) RELYM,RELYC
      READ(7,*) 
      READ(7,*) KSW
      READ(7,*) 
      READ(7,*) MAXIFB
      READ(7,*) ERRFB
      READ(7,*) 
      READ(7,*) RO,VISCOS
      READ(7,*)

      READ(7,*) EPS
      READ(7,*) 
      READ(7,*) MAXCPL
      READ(7,*) TNOEPS
      CLOSE(7)


      XLEN=RADCRY+RADCRU


      NNODX=2*NELEXC+1 + 2*NELEXM
      NNODXC=2*NELEXC+1
      NELEMM=(NELEXC+NELEXM)*NELEYM
      NNODM=NNODX*(2*NELEYM+1)-NELEMM


      OPEN (UNIT=8,FILE='initmesh.dat',STATUS='OLD',ERR=990)

      READ(8,*) NNODES,NELEMS
      CALL LENTST (NNODES,MXNOD,'MXNOD','MAIN')
      CALL LENTST (NELEMS,MXELE,'MXELE','MAIN')

      DO 100 NELE=1,NELEMS
         READ(8,*) ITYP,NODEL,(ELTOP(2+J,NELE),J=1,NODEL)
         ELTOP(1,NELE) = ITYP
         ELTOP(2,NELE) = NODEL
 100  CONTINUE

      DO 110 NODE=1,NNODES
         READ (8,*) (COORD(J,NODE),J=1,2)
 110  CONTINUE

      CLOSE(8)
 

      IPAR=NELEXC+NELEXM+1

      IBANDV=12*IPAR+7
      IDIAGV=6*IPAR+5
      LOV=6*IPAR+4
      IUPV=6*IPAR+4


      IBANDT=6*IPAR+3+2
      IDIAGT=3*IPAR+2+1
      LOT=3*IPAR+2
      IUPT=3*IPAR+2


c      WRITE(6,*) 'Initial velocity profile from file velocity.dat'
c      WRITE(6,*) '1 = yes, 0 = no'
c      READ(5,*) IFILE

      IFILE=0

      IF (IFILE .EQ. 1) THEN
         OPEN (UNIT=18,FILE='velocity.dat',STATUS='OLD',ERR=990)
         READ(18,*) NNODES
         DO 10 I=1,NNODES
            READ(18,601) NNN,XXX,YYY,UPREV(2*I-1),UPREV(2*I)
 10      CONTINUE
         CLOSE(18)

         OPEN (UNIT=18,FILE='temp.dat',STATUS='OLD',ERR=990)
         READ(18,*) NNODES
         DO 15 I=1,NNODES
            READ(18,*) TPREV(I)
 15      CONTINUE
         CLOSE(18)

      ELSE
         DO 20 I=1,NNODES
            TPREV(I)=TGAMIN
            UPREV(2*I-1)=0.0D0
            UPREV(2*I)=YVEL
 20      CONTINUE
      ENDIF

****************************************************************************
      ALIQ=TMELT
****************************************************************************


      ITER=1
      TIME=0.0D0

 6000 CONTINUE

      IND=NNODM-NNODX

      DO 120 I=1,NNODX
         SHCRU(I)=COORD(2,I)
         SHLIQ(I)=COORD(2,IND+I)
 120  CONTINUE
      CLOSE(12)

      DO 125 I=1,NNODXC
         FREEB(I)=SHLIQ(I)
 125  CONTINUE

      CALL VFILL (2*MXNOD,U,ZERO)
      CALL VFILL (MXNOD,T,ZERO)


      ITCOUP=1	

      WRITE(6,*) 
     +'Iteration,        time               rel.err. T      Rel.err. U'  

 7000 CONTINUE

      CALL LTIMST(DT,COORD,DIM,UPREV,NELEMS,NNODES,ELTOP,IELTOP,
     +     'DT','MAIN')

      TIME=TIME+DT

********* Calculation of temperatures **********************************

      CALL VFILL (2*MXNOD*IBMAX,A,ZERO)
      CALL VFILL (2*MXNOD,F,ZERO)

      DO 130 NELE=1,NELEMS
         NODEL = ELTOP(2,NELE)
         CALL ECOORD (NELE,GEOM,NODEL,ELTOP,IELTOP,COORD,ICOORD)
         CALL DIRECT (NELE,ELTOP,IELTOP,STEER)
         CALL CONVEL (NODEL,XVELO,YVELO,UPREV,STEER)

         DO 140 I=1,NODEL
            IP=STEER(I)
            XP=COORD(1,IP)
            YP=COORD(2,IP)

            XO=XP
            YO=YP
            
            IF (IP .LE. NNODM) THEN
               CALL MPOINT(XP,YP,XVELO(I),YVELO(I),DT,COORD,DIM,
     +              NNODX,NODEL,NODEL,ELEMCO,FLOWEL,IP,TPREV,
     +              IELTOP,ELTOP,XLEN,INDELE,SHCRU,SHLIQ,
     +              NNODM,NELEMM)
            ELSE
               CALL MPCRYS(XP,YP,XVELO(I),YVELO(I),DT,COORD,DIM,
     +              NODEL,ELEMCO,FLOWEL,IP,TPREV,
     +              IELTOP,ELTOP,INDELE,NNODM,NELEXC,NELE,I)
            ENDIF
            
            TELEM(I)=CONVT(XP,YP,XO,YO,NODEL,ELEMCO,FLOWEL,
     +           INDELE)
            
 140     CONTINUE
         
         CALL L2DFCV (ALOK,MLOK,FLOK,GEOM,NODEL,KSW,RO,CP,
     +        DK,TELEM,DT)
         
         CALL ASMBTE (A,ALOK,MLOK,F,FLOK,2*MXNOD,STEER,NODEL,
     +        IDIAGT)
         
 130  CONTINUE


      DO 150 NELE=NELEMM-NELEXM+1,NELEMM
         NODEL = ELTOP(2,NELE)
         SIDNUM=3
         CALL ECOORD (NELE,GEOM,NODEL,ELTOP,IELTOP,COORD,ICOORD)
         CALL DIRECT (NELE,ELTOP,IELTOP,STEER)
         CALL C2LELL (ALOK,FLOK,SIDNUM,GEOM,NODEL,KSW,ALFA,
     +        AREFT,DT)
         CALL VFILL(4*NLOK*NLOK,MLOK,ZERO)
         CALL ASMBTE (A,ALOK,MLOK,F,FLOK,2*MXNOD,STEER,NODEL,
     +        IDIAGT)
 150  CONTINUE
      
      DO 160 NELE=NELEMM+NELEXC,NELEMS,NELEXC
         NODEL = ELTOP(2,NELE)
         SIDNUM=2
         CALL ECOORD (NELE,GEOM,NODEL,ELTOP,IELTOP,COORD,ICOORD)
         CALL DIRECT (NELE,ELTOP,IELTOP,STEER)
         CALL C2LELL (ALOK,FLOK,SIDNUM,GEOM,NODEL,KSW,ALFA,AREFT,DT)
         CALL VFILL(4*NLOK*NLOK,MLOK,ZERO)
         CALL ASMBTE (A,ALOK,MLOK,F,FLOK,2*MXNOD,STEER,NODEL,
     +                IDIAGT)
 160  CONTINUE


      I=1
      DO 170 NELE=(NELEMM-(NELEXC+NELEXM)+1),(NELEMM-NELEXM)
         NODEL = ELTOP(2,NELE)
         SIDNUM=3
         DSIG=(FREEB(I+2)-FREEB(I))/(COORD(1,I+2)-COORD(1,I))
         ASTEF=1.0D0
         UEXT=RO*DLHEAT*YVEL/DSQRT(1.0D0+DSIG*DSIG)
         CALL ECOORD (NELE,GEOM,NODEL,ELTOP,IELTOP,COORD,ICOORD)
         CALL DIRECT (NELE,ELTOP,IELTOP,STEER)
         CALL C2LELL (ALOK,FLOK,SIDNUM,GEOM,NODEL,KSW,ASTEF,UEXT,DT)
         CALL VFILL(4*NLOK*NLOK,ALOK,ZERO)
         CALL VFILL(4*NLOK*NLOK,MLOK,ZERO)
         CALL ASMBTE (A,ALOK,MLOK,F,FLOK,2*MXNOD,STEER,NODEL,
     +                IDIAGT)
         I=I+2
 170  CONTINUE

      CALL BCONDT(A,F,2*MXNOD,NNODX,NNODXC,NNODES,NNODM,
     +     TGAMIN,TGAMOU,IDIAGT)


      CALL SOLVE(A,T,F,MXNOD,NNODES,LOT,IDIAGT,IUPT)



      CALL NEWIT(TPREV,T,NNODES,TERR)


***************  Calculation of velocities  ***********************************

         CALL VFILL (2*MXNOD*IBMAX,A,ZERO)
         CALL VFILL (2*MXNOD,F,ZERO)

         DO 200 NELE=1,NELEMM
            NODEL = ELTOP(2,NELE)
            CALL ECOORD (NELE,GEOM,NODEL,ELTOP,IELTOP,COORD,ICOORD)
            CALL DIRECT (NELE,ELTOP,IELTOP,STEER)
            DO 230 I=1,NODEL
               IP=STEER(I)
               XP=COORD(1,IP)
               YP=COORD(2,IP)

               XO=XP
               YO=YP

               CALL MPOINT(XP,YP,UPREV(2*IP-1),UPREV(2*IP),DT,COORD,
     +              DIM,NNODX,NODEL,2*NODEL,ELEMCO,FLOWEL,IP,UPREV,
     +              IELTOP,ELTOP,XLEN,INDELE,SHCRU,SHLIQ,
     +              NNODM,NELEMM)

               DO 240 J=1,NODEL
                  FLX(J)=FLOWEL(2*J-1)
                  FLY(J)=FLOWEL(2*J)
 240           CONTINUE
 
               UELEM(2*I-1)=CONVT(XP,YP,XO,YO,NODEL,ELEMCO,FLX,
     +              INDELE)
               UELEM(2*I)=CONVT(XP,YP,XO,YO,NODEL,ELEMCO,FLY,
     +              INDELE)

	       TELEM(I)=TPREV(ELTOP(2+I,NELE))
 230        CONTINUE

            CALL L2NAVI (ALOK,FLOK,GEOM,VISCOS,RO,DT,
     +           UELEM,EPS,NODEL,TELEM,GRAVX,GRAVY,HEXPCO,REFTEM)
            CALL ASSEMB(A,ALOK,F,FLOK,MXNOD,STEER,NODEL,IDIAGV)
 200     CONTINUE

         CALL BCOND(A,F,MXNOD,NNODX,NNODM,IDIAGV,COORD,ICOORD,
     +        XLEN,RADCRY,YVEL)
         CALL SOLVE(A,U,F,MXNOD,2*NNODM,LOV,IDIAGV,IUPV)

         DO 210 I=NNODM+1,NNODES
            U(2*I-1)=0.0D0
            U(2*I)=YVEL
 210     CONTINUE



         CALL NEWIT(UPREV,U,2*NNODES,VERR)

*************** End of calculation of velocities  ****************************

         WRITE(6,17) ITCOUP,TIME,TERR,VERR
 17      FORMAT(I9,3F20.8)

      IF ((VERR .GT. TNOEPS) .AND. (ITCOUP .LT. MAXCPL)) THEN
         ITCOUP=ITCOUP+1
         GOTO 7000
      ENDIF


      ILIN1=1
      ILIN2=1
      INDMEL=3*(NELEXC+NELEXM+1)-1
      IND1=0
      IND2=INDMEL
      DO 300 I=1,NNODXC
         ILIN1=1
         ILIN2=1
 400     IF ((T(I+IND1) .GT. TMELT) .AND. 
     +        (T(I+IND2) .GT. TMELT)) THEN


            IF ( ((I+IND2)+INDMEL) .GT. NNODM) THEN
               IF (ILIN2 .EQ.1) THEN
                  IND2=IND2+NNODX+(NELEXC+1)
                  ILIN2=2
               ELSE
                  IND2=IND2+3*(NELEXC+1)-1
               ENDIF
            ELSE
               IND2=IND2+INDMEL
            ENDIF

            IF ( ((I+IND1)+INDMEL) .GT. NNODM) THEN
               IF (ILIN1 .EQ.1) THEN
                  IND1=IND1+NNODX+(NELEXC+1)
                  ILIN1=2
               ELSE
                  IND1=IND1+3*(NELEXC+1)-1
               ENDIF
            ELSE
               IND1=IND1+INDMEL
            ENDIF

            GOTO 400
         ELSE
            FREEBO(I)=FREEB(I)
            FREEB(I)=(TMELT-T(I+IND2))*(COORD(2,I+IND1)-
     +           COORD(2,I+IND2))/
     +           (T(I+IND1)-T(I+IND2))+COORD(2,I+IND2)
         ENDIF
 300  CONTINUE

      write(6,*) '*** Free boundary iteration'
      write(6,411) freeb(1),freeb(3),freeb(5),freeb(7),freeb(9),
     +     freeb(11),freeb(13),ENORM(FREEB,FREEBO,NNODXC)

 411  format(8f9.5)

      OPEN(11,FILE='velocity.dat',STATUS='UNKNOWN',ERR=991)
      write(11,*) nnodes
      do 610 i=1,nnodes
         write(11,601) i,coord(1,i),coord(2,i),u(2*i-1),u(2*i)
 610  continue
      write(11,*) nelexc,nelexm
      write(11,*) neleyc,neleym
 601  format(i5,5f14.9)
      CLOSE(11)

      OPEN(UNIT=12,FILE='temp.dat',STATUS='UNKNOWN')
      write(12,*) nnodes
      do 529 i=1,nnodes
         write(12,*) t(i)
 529  continue
      CLOSE(12)

      OPEN (UNIT=8,FILE='mesh.dat',STATUS='UNKNOWN',ERR=991)
      write(8,*) nnodes,nelems
      do 509 nele=1,nelems
         write(8,*) (ELTOP(J,NELE),J=1,NODEL+2)
 509  continue
      do 519 i=1,nnodes
         write(8,*) (COORD(J,i),J=1,2)
 519  continue
      CLOSE(8)


      IF ((ITER .LT. MAXIFB) .AND. 
     +     (ENORM(FREEB,FREEBO,NNODXC) .GT. ERRFB)) THEN

         CALL VFILL(DIM*MXNOD,COORD,ZERO)
 
         CALL MESHG(FREEB,COORD,ICOORD,NELEXC,NELEXM,NELEYC,NELEYM,
     +        XLEN,YLEN,RELXC,RELXM,RELYM,RELYC,RO,GASDEN,SURTEN,COANG,
     +        RADCRY,RADCRU)

         ITER=ITER+1
         GOTO 6000
      ENDIF

      GOTO 1000
 990  CONTINUE
      print *,'Error in input files'
      STOP
 991  CONTINUE
      print *,'Error in output files' 
      STOP
 1000 CONTINUE

      END                                 
