I tend to use temporary variables to “ignore” some subroutine outputs. Since Fortran doesn’t have a command like tilde (~) in Matlab, I have to get the outputs fron the subroutine but I assign them to the same temporary variable with the right size. It has been my preference to make things look cleaner, no practical purposes.

So far for years I had no problems with it; today, I am getting NaN. I isolated the problem to a matrix inverse operation in the “r_calc” subroutine but I am not sure what is happening. If I replace that with a transpose (which actually gives the same result for this matrix), the problem disappears.

My question is, is it bad practive to assign multiple subroutine outputs to the same variable in the main code? Why am I having an issue with that when I use matrix inverse? I appreciate any help.

Below is a minimum working example:

program example

implicit none

```
real(kind=8) :: temp1(3,3)
real(kind=8) :: temp2(3,3)
real(kind=8) :: temp3(3,3)
real(kind=8) :: temp4(3,3)
real(kind=8) :: temp5(3,3)
real(kind=8) :: r_PN_U(3)
call r_calc(4.2d0, &
0d0, 0d0, &
0d0, 0d0, &
0d0, 0d0, &
0d0, &
0d0, 0d0, &
0d0, &
0d0, 0d0, &
0d0, 0d0, &
r_PN_U, &
temp1, temp1, temp1, temp1, temp1)
print *, r_PN_U ! gives NaN in mat_inv_3 is used in r_calc subroutine
call r_calc(4.2d0, &
0d0, 0d0, &
0d0, 0d0, &
0d0, 0d0, &
0d0, &
0d0, 0d0, &
0d0, &
0d0, 0d0, &
0d0, 0d0, &
r_PN_U, &
temp1, temp2, temp3, temp4, temp5)
print *, r_PN_U ! gives correct values even with mat_inv_3
```

end program example

subroutine r_calc(rin, &

e, ed, &

f, fd, &

v, vd, &

vp, &

w, wd, &

wp, &

theta, thetad, &

y_pos, z_pos, &

r_PN_U, &

C_DU_theta, C_DU_beta, C_DU_zeta, C_DU, C_UD)

implicit none

```
real(kind=8), intent(in) :: rin
real(kind=8), intent(in) :: e, ed
real(kind=8), intent(in) :: f, fd
real(kind=8), intent(in) :: v, vd
real(kind=8), intent(in) :: vp
real(kind=8), intent(in) :: w, wd
real(kind=8), intent(in) :: wp
real(kind=8), intent(in) :: theta, thetad
real(kind=8), intent(in) :: y_pos, z_pos
real(kind=8), intent(out) :: r_PN_U(3)
real(kind=8), intent(out) :: C_DU_theta(3,3)
real(kind=8), intent(out) :: C_DU_beta (3,3)
real(kind=8), intent(out) :: C_DU_zeta (3,3)
real(kind=8), intent(out) :: C_DU (3,3)
real(kind=8), intent(out) :: C_UD (3,3)
real(kind=8) :: beta, zeta
beta = -atan(wp) ! [rad], flap down angle
zeta = atan(vp) ! [rad], lead angle
call angle2dcm(theta, beta, zeta, C_DU_theta, C_DU_beta, C_DU_zeta, C_DU)
call mat_inv_3(C_DU, C_UD) ! results in NaN in r_PN_U output
! C_UD = transpose(C_DU) ! gives the same result as inverse, eliminates the NaN issue
r_PN_U = [rin+e+f, v, w] + matmul(C_UD, [0d0, y_pos, z_pos])
```

end subroutine r_calc

subroutine angle2dcm(phi, theta, psi, C_phi, C_theta, C_psi, C_out)

implicit none

```
! Calculates the direction cosine matrix in psi - theta - phi (3 - 2 - 1) order
! Difference from "angle2dcm" subroutine is the extra outputs
real(kind=8), intent(in) :: phi, theta, psi
real(kind=8), intent(out) :: C_psi(3,3), C_theta(3,3), C_phi(3,3), C_out(3,3)
C_phi(1,1:3) = [1d0, 0d0, 0d0]
C_phi(2,1:3) = [0d0, cos(phi), sin(phi)]
C_phi(3,1:3) = [0d0, -sin(phi), cos(phi)]
C_theta(1,1:3) = [cos(theta), 0d0, -sin(theta)]
C_theta(2,1:3) = [ 0d0, 1d0, 0d0]
C_theta(3,1:3) = [sin(theta), 0d0, cos(theta)]
C_psi(1,1:3) = [ cos(psi), sin(psi), 0d0]
C_psi(2,1:3) = [-sin(psi), cos(psi), 0d0]
C_psi(3,1:3) = [ 0d0, 0d0, 1d0]
C_out = matmul(C_phi, matmul(C_theta,C_psi)) ! psi - theta - phi (3 - 2 - 1) order
```

end subroutine angle2dcm

subroutine mat_inv_3(A, B)

implicit none

```
real(kind=8), intent(in) :: A(3,3)
real(kind=8), intent(out) :: B(3,3)
real(kind=8) :: det
det = 1d0/(A(1,1)*A(2,2)*A(3,3) - A(1,1)*A(2,3)*A(3,2)&
- A(1,2)*A(2,1)*A(3,3) + A(1,2)*A(2,3)*A(3,1)&
+ A(1,3)*A(2,1)*A(3,2) - A(1,3)*A(2,2)*A(3,1))
B(1,1) = +det * (A(2,2)*A(3,3) - A(2,3)*A(3,2))
B(2,1) = -det * (A(2,1)*A(3,3) - A(2,3)*A(3,1))
B(3,1) = +det * (A(2,1)*A(3,2) - A(2,2)*A(3,1))
B(1,2) = -det * (A(1,2)*A(3,3) - A(1,3)*A(3,2))
B(2,2) = +det * (A(1,1)*A(3,3) - A(1,3)*A(3,1))
B(3,2) = -det * (A(1,1)*A(3,2) - A(1,2)*A(3,1))
B(1,3) = +det * (A(1,2)*A(2,3) - A(1,3)*A(2,2))
B(2,3) = -det * (A(1,1)*A(2,3) - A(1,3)*A(2,1))
B(3,3) = +det * (A(1,1)*A(2,2) - A(1,2)*A(2,1))
```

end subroutine