请帮忙看一下这个子程序

最近在做一个坐标变换后焊接热源的子程序，焊接路径不是单纯沿着坐标轴方向运动。大家帮忙看一下下面编译的子程序。现在比较迷惑的是子程序里面的焊接路径这么编写对不对，先谢过了！
subroutine flux(f,temflu,mibody,time)
implicit real*8 (a-h,o-z)
dimension f(2),mibody(*),temflu(*)
c* * * * * *
c
c user subroutine for non-uniform flux input.
c
c f(1) flux value (to be defined)
c f(2) derivative of flux with respect to temperature
c (to be defined; optional, this might improve the
c convergence behavior)
c
c temflu(1) estimated temperature
c temflu(2) previous volumetric flux
c temflu(3) temperature at beginning of increment
c temflu(4,5,6)integration point coordinates
c mibody(1) element number
c mibody(2) flux type
c mibody(3) integration point number
c mibody(4) flux index
c mibody(5) not used
c mibody(6) =1 : heat transfer
c =2 : joule
c =3 : bearing
c =4 : electrostatic
c =5 : magnetostatic
c =6 : acoustic
c mibody(8) layer number for heat transfer shells elements
c and volume flux
c time time
c
c* * * * * *
f(1)=0 !定义初始热流
a=0.0055 !x半轴长
b=0.008 !y半轴长
c1=0.013 !z前半轴长
c2=0.013 !z后半轴长
r=0.002 !旋转半径
vx=0.00383 !x方向速度
vz=0.1 !z方向速度
x=temflu(4)
y=temflu(5)
z=temflu(6)
x0=vx*time !定义焊接路径
z0=vz*time
x1=temflu(4)-x0 !x变换
y1=temflu(5)
z1=temflu(6)-z0 !z变换
if(z1.ge.0) then
f=2e6*exp(-3*((x1/a)**2+(y1/b)**2+(z1/c1)**2))
else if(z1.lt.0) then
f=3e6*exp(-3*((x1/a)**2+(y1/b)**2+(z1/c2)**2))
else
f=0
end if
return
end