SUBROUTINE DFLUX(FLUX,SOL,KSTEP,KINC,TIME,NOEL,NPT,COORDS,
     1 JLTYP,TEMP,PRESS,SNAME)
C
      INCLUDE 'ABA_PARAM.INC'
C
      DIMENSION FLUX(2), TIME(2), COORDS(3)
      CHARACTER*80 SNAME
!      
!----------------------------------------------------------------------         
!
!      Declaration of variables       
       REAL :: f1, f2
       REAL :: ALP, Rf, T1, Imax, tmax, Laser_I
       REAL :: coordx, coordy, coordy2, abldepth
       PARAMETER (OriSurfy = 0.00999999978)
       PARAMETER (nsurfeles = 94)
       PARAMETER (nsurfnode = 95)
       REAL :: Val_abldepth,r,b,qflux,x,x1
       PARAMETER (surfeleini = 7897)
       PARAMETER (surfeleinc = 1)
       PARAMETER (nlayer = 84)
       PARAMETER (layerinc = 94)
       PARAMETER (beta = 1.8E-01) 
       PARAMETER (atoms = 4.4803895E-026)
       PARAMETER (boltz = 1.38064852E-023)
       PARAMETER (tboil = 2743)
       PARAMETER (pboil = 0.1E06)
       PARAMETER (lvheat = 10.78E06)
!      PARAMETER (lmelt = 3.96E05)
       PARAMETER (critemp = 8860)
       PARAMETER (meltemp = 933)
       PARAMETER (critrho = 300)
       
       REAL :: rho,n,recesrate,rho_2
       REAL :: coef1, coef2, Talp
!
!----------------------------------------------------------------------  
!     
      REAL, DIMENSION (nsurfnode) :: XCOORD, YCOORD
      COMMON /COORDData/ XCOORD, YCOORD  
!      
!----------------------------------------------------------------------               

       IF (JLTYP == 0) THEN
          DO I=1,nsurfeles
          J=surfeleini+surfeleinc*(I-1)
          IF (NOEL == J .AND. NPT == 2) THEN
              XCOORD(I) = COORDS(1)
              XCOORD(I) = XCOORD(I) *1.0E00
              YCOORD(I) = -COORDS(2) + OriSurfy
              YCOORD(I) = YCOORD(I) *1.0E00
          ENDIF
          END DO
          k=surfeleini+surfeleinc*(nsurfeles-1)
          IF (NOEL == k .AND. NPT == 1) THEN
              XCOORD(nsurfnode) = COORDS(1)
              XCOORD(nsurfnode) = XCOORD(nsurfnode) *1.0E00
              YCOORD(nsurfnode) = -COORDS(2) + OriSurfy
              YCOORD(nsurfnode) = YCOORD(nsurfnode) *1.0E00
          ENDIF
          
      n = 0.333333333      
          
      IF (SOL <= meltemp) THEN
          rho = 2700
      ELSEIF (SOL > meltemp .and. SOL <= critemp) THEN
          rho = critrho*(1+0.75*(1-SOL/critemp)
     $     +3*(1-SOL/critemp)**n)
      ELSEIF (SOL > critemp) THEN
          rho = critrho
      ENDIF         
          
      recesrate = (1-beta)*sqrt(atoms/(6.28*boltz*SOL))*pboil/rho*
     $   exp((atoms*lvheat)/boltz*(1/tboil-1/SOL))
      recesrate = recesrate * 1.0E3
      
      rho_2 = rh0 * 1.0E-9
          
      IF (SOL < tboil) THEN
          FLUX(1) = 0
      ELSEIF (SOL >= tboil) THEN
          FLUX(1) = -lvheat*rho_2*recesrate
      ENDIF
           
       ELSEIF (JLTYP == 1) THEN
!
!----------------------------------------------------------------------          
!
!     Get the orginal y coordinates      
!      
      coordx = COORDS(1)
      coordy = -COORDS(2) + OriSurfy
! 
!---------------------------------------------------------------------- 
      
       DO I=1,nsurfeles
          J=surfeleini+surfeleinc*(I-1)
          DO K=0,nlayer
          IF (NOEL == J-layerinc*K .AND. NPT == 3) THEN
              coordy2 = coordy - YCOORD(I)
          ELSEIF (NOEL == J-layerinc*K .AND. NPT == 1) THEN
              coordy2 = coordy - YCOORD(I)
          ELSEIF (NOEL == J-layerinc*K .AND. NPT == 4) THEN
              coordy2 = coordy - YCOORD(I+1)     
          ELSEIF (NOEL == J-layerinc*K .AND. NPT == 2) THEN
              coordy2 = coordy - YCOORD(I+1)                 
          ENDIF
          END DO
       END DO
       
      if (noel == 7990) then
      print *,kinc,npt,coordy,coordy2
      endif    
   
! 
!----------------------------------------------------------------------      
!
!         Parameters
!      
!      coef1 = 4.0846153846E-02
!      coef2 = 3.6003615385E+02
      Talp = critemp
      
      IF (SOL <= meltemp) THEN 
          ALP = 3.83E05
      ELSEIF (SOL > meltemp .and. SOL <= 6031) THEN
          ALP = -2.6733633367E+04 * SOL + 4.0825265685E+08   
          ALP = ALP * 1E-03 
      ELSEIF (SOL > 6031 .and. SOL <= Talp) THEN
          ALP = -8.5962888740E+04 * SOL + 7.6546429600E+08
          ALP = ALP * 1E-03
      ELSE
          ALP = 3.83E03
      ENDIF
      
!      Rf = 4.6E-01
      IF (SOL <= Talp) THEN
          Rf = -4.4305446964E-05 * SOL + 8.1329163409E-01
      ELSE
          RF = 0.46
      ENDIF
      
           
      T1 = TIME(1)
      Imax = 3.78E07
      tmax = 2.5E-09 
      Laser_I = Imax*((T1/tmax)**7)*EXP(7*(1-T1/tmax))
      r=50E-3
      b=-LOG(0.1)/((0.55*r)**2)     
      x1 = 75.E-03
      x = COORDS(1)      
!      
!---------------------------------------------------------------------- 
!
      
      IF (coords(1) > 25.0E-03) THEN
          f1 = (1-Rf)*ALP*EXP(-ALP*coordy2)
          Qflux = f1*Laser_I
          shape=exp(-b*(x-x1)**2)
          FLUX(1) = Qflux*shape
      ENDIF
      
      ENDIF
      
!      
      RETURN
      END SUBROUTINE DFLUX
!      
!----------------------------------------------------------------------
!----------------------------------------------------------------------
!     
      SUBROUTINE umeshmotion(uref,ulocal,node,nndof,lnodetype,alocal,
     $     ndim,time,dtime,pnewdt,kstep,kinc,
     $     kmeshsweep,jmatyp,jgvblock,lsmooth)
!
      include 'ABA_PARAM.INC'
!
!     umeshmotion subroutine
      dimension ulocal(ndim)
      dimension alocal(ndim,*),time(2)
      dimension jmatyp(*),jgvblock(*)

!     getnodetoelemconn routine
      parameter (maxnelems = 10000)
      dimension jelemlist(maxnelems),jelemtype(maxnelems)

!     getvrmavgatnode routine
      dimension array(15)

!     getpartinfo routine
      character*80 partname
      character*80 cpname
!      
!----------------------------------------------------------------------  
!
      REAL :: x1, x2, n
      REAL :: tmp1, tmp2
      REAL :: xxx, xx1, crd1, crd2
      REAL :: rho
!
!----------------------------------------------------------------------
!     
!     coefficient for hertz-knudeson equation
      PARAMETER (beta = 1.8E-01) 
      PARAMETER (atoms = 4.4803895E-026)
      PARAMETER (boltz = 1.38064852E-023)
      PARAMETER (tboil = 2743)
      PARAMETER (pboil = 0.1E06)
      PARAMETER (lvheat = 10.78E06)
      PARAMETER (critemp = 8860)
      PARAMETER (meltemp = 933)
      PARAMETER (critrho = 300)
      PARAMETER (nsurfnodes = 95)
      PARAMETER (OriSurfy = 0.00999999978)      
!      
!---------------------------------------------------------------------- 
!       
      LOCNUM = 0
      JRCD = 0
      PARTNAME = ' '
      CHARLENGTH = UREF
      JTYP = 1
      
      ndt1 = NODE
      ndx1 = NODE
      ndt2 = ndt1 - 1
      ndx2 = ndx1 - 1
      ltrn = 0
c
      call GETVRN(ndt1,'NT',array,jrcd,JGVBLOCK,ltrn)
      tmp1 = array(1)
      call GETVRN(ndx1,'COORD',array,jrcd,JGVBLOCK,ltrn)
      crd1 = array(2)
      call GETVRN(ndt2,'NT',array,jrcd,JGVBLOCK,ltrn)
      tmp2 = array(1)
      call GETVRN(ndx2,'COORD',array,jrcd,JGVBLOCK,ltrn)
      crd2 = array(2)
      
!      
!---------------------------------------------------------------------- 
!            
      n = 0.333333333
      
!----------------------------------------------------------------------            
!      
      if (tmp1 <= meltemp) then
          rho = 2700
      elseif (tmp1 > meltemp .and. tmp1 <= critemp) then
          rho = critrho*(1+0.75*(1-tmp1/critemp)
     $     +3*(1-tmp1/critemp)**n)
      elseif (tmp1 > critemp) then
          rho = critrho
      endif
!      
!----------------------------------------------------------------------            
!      
          x1 = 0
!          x2 = (1-beta)*sqrt(atoms/(6.28*boltz*tmp1))*pboil/rho* 
!     $         exp(lvheat/boltz*(1/tboil-1/tmp1))          
          x2 = (1-beta)*sqrt(atoms/(6.28*boltz*tmp1))*pboil/rho*
     $         exp((atoms*lvheat)/boltz*(1/tboil-1/tmp1))
          
!         unit conversion from m to mm          
          x2 = x2 * 1E03
!
          ULOCAL(1) = ULOCAL(1) - ALOCAL(1,1)*x1 - ALOCAL(2,1)*x2
          ULOCAL(2) = ULOCAL(2) - ALOCAL(1,2)*x1 - ALOCAL(2,2)*x2
                   
!      
!----------------------------------------------------------------------          
      
      LSMOOTH = 1
      
!
      RETURN
      END SUBROUTINE umeshmotion    
!----------------------------------------------------------------------