      SUBROUTINE MESHG(FREEB,COORD,ICOORD,NELEXC,NELEXM,NELEYC,NELEYM,
     +     XLEN,YLEN,RELXC,RELXM,RELYM,RELYC,RO,ROGAS,SIGMA,ALFA,
     +     RADCRY,RADCRU)
      INTEGER ICOORD,NELEXC,NELEXM,NELEYC,NELEYM
      DOUBLE PRECISION FREEB(*),COORD(ICOORD,*),XLEN,YLEN,
     +     RELXC,RELXM,RELYM,RELYC,RO,ROGAS,SIGMA,ALFA,RADCRY,RADCRU

      INTEGER NNODX
      INTEGER NPOINT,NXMAX,KK,NY,J
      PARAMETER(NXMAX=100)
      DOUBLE PRECISION XCOORD(NXMAX),YCOORD(NXMAX),YTOP(NXMAX),
     +     YBOT(NXMAX),YBOUND(NXMAX),YAPU(NXMAX),B(NXMAX),YMID(NXMAX)
      DOUBLE PRECISION RODIF,G,BETA,DIMEN,YSTEP,Y,X,YTPP,YSCL,SINIT,
     +     THETA,CX,CX0,A,XORIG,RR,X0

      PARAMETER(THETA=1.0D-6)

      DOUBLE PRECISION LIQINF,densit

************************************
      densit=2.49d3
************************************



      NNODX=NELEXC+NELEXM+1

      DO 10 I=1,NELEXC+1
        YBOUND(I)=FREEB(2*I-1)
        YBOT(I)=0.0D0
        YTOP(I)=YLEN
 10   CONTINUE


      SINIT=YBOUND(NELEXC+1)
      RODIF=densit-ROGAS
      G=9.82D0
      BETA=2.0D4*SIGMA/(RODIF*G)
      DIMEN=DSQRT(2.0D0/BETA)


         RR=RADCRY*DIMEN*1.0D2

         Y=1.0D-18
         X0=1.0D12
         YSTEP=0.1d0
         NY=0
         NPOINT=100

         J=1
         DO 220 I=1,NPOINT
            X=LIQINF(ALFA,RR,BETA,Y,NY)
            IF ((NY .EQ. 1) .OR. (X .GT. X0) .OR. (X .LT. RR)) GOTO 221
            XCOORD(J)=X/(DIMEN*1.0D2)
            YCOORD(J)=Y/(DIMEN*1.0D2)
            Y=Y+YSTEP
            J=J+1
            X0=X
 220     CONTINUE
 221     CONTINUE
         NPOINT=J-1

         IF (XCOORD(1) .LT. RADCRY+RADCRU) THEN
            print *,'Length of meniscus insufficient'
            STOP
         ENDIF

         YTPP=YCOORD(NPOINT)+(YCOORD(NPOINT)-YCOORD(NPOINT-1))/
     +        (XCOORD(NPOINT)-XCOORD(NPOINT-1))*(RADCRY-XCOORD(NPOINT))

         YSCL=SINIT-YTPP
         XCOORD(NPOINT+1)=RADCRY
         YCOORD(NPOINT+1)=SINIT

         DO 225 I=1,NPOINT
            YCOORD(I)=YCOORD(I)+YSCL
 225     CONTINUE

      CX=RADCRY
      DO 20 I=NELEXC+2,NNODX
         IF (ABS(RELXM-1.0D0) .LT. THETA) THEN	
            CX=CX+RADCRU/DBLE(NELEXM)
         ELSE IF (I .EQ. NELEXC+2) THEN
            A=RADCRU*(1.0D0-RELXM)/(1.0D0-RELXM**DBLE(NELEXM))	
            CX=CX+A
         ELSE
            CX=CX+RELXM*A
            A=RELXM*A	
         ENDIF		
         
         DO 25 J=NPOINT+1,2,-1
            
            IF ((XCOORD(J) .LT. CX) .AND. 
     +           (CX .LT. XCOORD(J-1))) THEN
               YBOUND(I)=(YCOORD(J-1)-YCOORD(J))/
     +              (XCOORD(J-1)-XCOORD(J))*(CX-XCOORD(J-1))+
     +              YCOORD(J-1)
               GOTO 26
            ENDIF
            
 25      CONTINUE
 26      CONTINUE
         YBOT(I)=0.0D0
         YTOP(I)=YLEN
 20   CONTINUE
      
      XORIG=0.0D0
      CX=0.0D0
      KK=1
      
      DO 140 I=1,2*NELEYM+1


         DO 150 J=1,NELEXC+1
            IF (MOD(I-1,2) .EQ. 0) THEN
               COORD(1,KK)=CX
               COORD(2,KK)=YBOT(J)
            ELSE
               COORD(1,KK)=CX
               COORD(2,KK)=(YBOT(J)+YAPU(J))/2.0D0
            ENDIF

            CX0=CX
            YAPU(J)=YBOT(J)

            IF (J .LT. NELEXC+1) THEN
               IF (ABS(RELXC-1.0D0) .LT. THETA) THEN	
                  CX=CX+RADCRY/DBLE(NELEXC)
               ELSE IF (J .EQ. 1) THEN
                  A=RADCRY*(1.0D0-RELXC)/(1.0D0-RELXC**DBLE(NELEXC))
                  CX=CX+A
               ELSE
                  CX=CX+RELXC*A
                  A=RELXC*A	
               ENDIF		
            ENDIF

            IF (MOD(I-1,2) .EQ. 0) THEN
               IF (J .LT. NELEXC+1) THEN
                  KK=KK+1
                  COORD(1,KK)=(CX+CX0)/2.0D0
                  COORD(2,KK)=(YBOT(J)+YBOT(J+1))/2.0D0
               ENDIF
            ENDIF
            KK=KK+1

 150     CONTINUE

         DO 155 J=NELEXC+2,NNODX
            CX0=CX


	    IF (ABS(RELXM-1.0D0) .LT. THETA) THEN	
               CX=CX+RADCRU/DBLE(NELEXM)

	    ELSE IF (J .EQ. NELEXC+2) THEN
	       A=RADCRU*(1.0D0-RELXM)/(1.0D0-RELXM**DBLE(NELEXM))	
	       CX=CX+A
	    ELSE
	       CX=CX+RELXM*A
	       A=RELXM*A	
	    ENDIF		

            IF (MOD(I-1,2) .EQ. 0) THEN
               COORD(1,KK)=(CX+CX0)/2.0D0
               COORD(2,KK)=(YBOT(J-1)+YBOT(J))/2.0D0
               KK=KK+1
            ENDIF

            IF (MOD(I-1,2) .EQ. 0) THEN
               COORD(1,KK)=CX
               COORD(2,KK)=YBOT(J)
            ELSE
               COORD(1,KK)=CX
               COORD(2,KK)=(YBOT(J)+YAPU(J))/2.0D0
            ENDIF

            KK=KK+1

            YAPU(J)=YBOT(J)
 155     CONTINUE


         CX=XORIG

         IF (MOD(I-1,2) .EQ. 0) THEN
            DO 160 K=1,NNODX
               IF (ABS(RELYM-1.0D0) .LT. THETA) THEN	
                  YBOT(K)=YBOT(K)+YBOUND(K)/DBLE(NELEYM)
               ELSE IF (I .EQ. 1) THEN
                  B(K)=YBOUND(K)*(1.0D0-RELYM)/
     +                 (1.0D0-RELYM**DBLE(NELEYM))	
                  YBOT(K)=YBOT(K)+B(K)
               ELSE
                  YBOT(K)=YBOT(K)+RELYM*B(K)
                  B(K)=RELYM*B(K)	
               ENDIF		

 160        CONTINUE   
         ENDIF
 140  CONTINUE


         DO 165 I=1,NELEXC+1
            YMID(I)=YBOUND(I)
 165     CONTINUE

      DO 170 I=1,2*NELEYC
 
         IF (MOD(I-1,2) .EQ. 0) THEN
            DO 180 K=1,NELEXC+1

               IF (ABS(RELYC-1.0D0) .LT. THETA) THEN	
                  YAPU(K)=YAPU(K)+(YTOP(K)-YBOUND(K))/DBLE(NELEYC)
               ELSE IF (I .EQ. 1) THEN
                  B(K)=(YTOP(K)-YBOUND(K))*(1.0D0-RELYC)/
     +                 (1.0D0-RELYC**DBLE(NELEYC))	
                  YAPU(K)=YAPU(K)+B(K)
               ELSE
                  YAPU(K)=YAPU(K)+RELYC*B(K)
                  B(K)=RELYC*B(K)	
               ENDIF		

 180        CONTINUE   
         ENDIF
         DO 190 J=1,NELEXC+1

            IF (MOD(I,2) .EQ. 0) THEN
               COORD(1,KK)=CX
               COORD(2,KK)=YAPU(J)
            ELSE
               COORD(1,KK)=CX
               COORD(2,KK)=(YAPU(J)+YMID(J))/2.0D0
            ENDIF

            KK=KK+1
            CX0=CX
            YMID(J)=YAPU(J)

            IF (ABS(RELXC-1.0D0) .LT. THETA) THEN	
               CX=CX+RADCRY/DBLE(NELEXC)
            ELSE IF (J .EQ. 1) THEN
               A=RADCRY*(1.0D0-RELXC)/(1.0D0-RELXC**DBLE(NELEXC))	
               CX=CX+A
            ELSE
               CX=CX+RELXC*A
               A=RELXC*A	
            ENDIF		
            IF ((MOD(I,2) .EQ. 0) .AND. (J .LT. NELEXC+1)) THEN
               COORD(1,KK)=(CX+CX0)/2.0D0
               COORD(2,KK)=(YAPU(J)+YAPU(J+1))/2.0D0
               KK=KK+1
            ENDIF

 190     CONTINUE
         CX=XORIG
 170  CONTINUE

      RETURN
      END
