[流程控制] 新人求帮助，热源模型参数循环加载问题

Program www_fcode_cn
    Implicit None
    Real(8) :: x(4,2),y
    Integer :: iDo,i
    Real(8), External   :: F
    !    ! 设置您的第一个参数
    !
    x(1,1) = ?
    x(2,1) = ?
    x(3,1) = ?
    x(4,1) = ?
    !
    ! 设置您的第二个参数
    !
    x(1,2) = ?
    x(2,2) = ?
    x(3,2) = ?
    x(4,2) = ?
    iDo = 0
    Do 
        iDo = iDo + 1
        Do i = 1 , 2
            y = F(X(:,i))
            write(*,*) '第',iDo,'次循环的第',i,'个参数计算结果为：',y
        EndDo
    EndDo  
End Program www_fcode_cn

FUNCTION F(X)
    Implicit None
    real(8), Intent(in) :: X(4)
    real(8) :: XX,YY,ZZ,TT
    real(8) :: QF,QR,AF,AR,B,C,X0,Y0,Z0,VY,AY
    Integer :: M1
    Real(8) :: PIDEG
    Real(8) :: XD,YD,ZD,SA,CA,A1,A2,XL,YL,ZL,F,V1,V2,V3,QC,AC,COND
!
!   F = QC * V1 * V2 * V3 with
!   V1 = exp( -( YY-Y0-VY*TT )^2/AC^2 )
!   V2 = exp( -( XX-X0 )^2/B^2 )
!   V3 = exp( -( ZZ-Z0 )^2/C^2 )
!   if ( -YY + Y0 +VY*TT ) greater than 0
!     QC = QF et AC = AF
!   else
!     QC = QR et AC = AR
!
!
! Input
!
      XX = X(1)   !; X Coordinate
      YY = X(2)   !; Y Coordinate
      ZZ = X(3)   !; Z Coordinate
      TT = X(4)   !; Time
!
! Variables
!
      QF = 72.1725997925    !; Maximal front source intensity
      QR = 36.0862998962    !; Maximal rear  source intensity
      AF = 3    !; Gaussian parameter
      AR = 4    !; Gaussian parameter
      B  = 2    !; Gaussian parameter
      C  = 3    !; Gaussian parameter
      X0 = 0    !; X initial location of source center
      Y0 = 0    !; Y initial location of source center
      Z0 = 0    ! Z initial location of source center
      VY = 0    !; Source displacement velocity
      AY = 0    !; Angle of torch [deg.]
!
! Constant
!
      M1 = -1
      PIDEG = ATAN(1.)
      PIDEG = PIDEG / 45.
      AY = AY * PIDEG
!
! Transformation of global to local coordinates
!
      XD   = XX - X0
      YD   = VY * TT
      YD   = YD + Y0
      ZD   = ZZ - Z0
!
! Source rotation about Y axis 
!
      SA   = SIN( AY )
      SA   = - SA
      CA   = COS( AY )
      A1   = XD * CA
      A2   = ZD * SA
      XL   = A1 + A2
      YL   = YY - YD
      A1   = ZD * CA
      A2   = XD * SA
      ZL   = A1 - A2
!
! Condition computation, QC and AC initialisation
!
      COND = VY * YL
      IF (VY .EQ. 0.) COND = YL
      QC   = QR
      AC   = AR
      IF( COND .GT. 0. ) QC = QF
      IF( COND .GT. 0. ) AC = AF
!
! V1 computation
!
      A1   = YL * YL
      A2   = AC * AC
      A2   = A1 / A2
      A2   = M1 * A2
      V1   = EXP( A2 )
!
! V2 computation
!
      A1   = XL * XL
      A2   = B * B
      A2   = A1 / A2
      A2   = M1 * A2
      V2   = EXP( A2 )
!
! V3 computation
!
      A1   = ZL * ZL
      A2   = C * C
      A2   = A1 / A2
      A2   = M1 * A2
      V3   = EXP( A2 )
!
! F computation
!
      F    = QC * V1
      F    = F * V2
      F    = F * V3
      RETURN
      END