Dear NREL,
I am looking to extract some data for AeroDyn for input into another solver. I have a couple of questions regarding axis system and conventions etc. which do not seem clearly explained in the code or documentation.
BladeLoad contains forces and moments on the rotating blade as a function of time.

Are these loads in the global AeroDyn axis system? If not, can you provide information (i.e. transformations) regarding the local axis system (e.g. origin and axes) which they relate to? I am using a simple model with straight blades and steady wind as a test case, and it looks like force(1), force(2) and moment(3) contain nonzero values while the other terms are all zero. Hence, I suspect they are local to the blade.

Also, I tried to lump the loads in point mesh (I guess similar to BladePtLoads but collocated with BladeLoad) and I observed the moments components are not zero which I am not sure why?

In the AeroDyn outputs calculation these forces and moments are also postmultiplied by an array called WithoutSweepPitchTwist. Does this represent another transformation out of the global system?
y%BladeLoad(k)%Force(:,j) = matmul( force, m%WithoutSweepPitchTwist(:,:,j,k) ) ! force per unit length of the jth node in the kth blade
y%BladeLoad(k)%Moment(:,j) = matmul( moment, m%WithoutSweepPitchTwist(:,:,j,k) ) ! moment per unit length of the jth node in the kth blade

Documentation states that these loads are per unit length. Can you provide some information about how these loads are integrated along the blade element length to produce point loads (e.g. which are suitable for application to nodes in ElastoDyn and analogous to the coupling I am trying)?

Also, in the definition of the BladeMotion and BladeRootMotion orientations, I am defining a DCM matrix based on the actual nodal rotations of my elements. This approach is different to the one used in the driver, for example for the BladeRootMotion:
theta(1) = (k1)*TwoPi/real(DvrData%numBlades,ReKi)
theta(2) = DvrData%precone
theta(3) = blade_pitch(k)
orientation = EulerConstruct(theta)
AD%u(1)%BladeRootMotion(k)%Orientation( :,:,1) = matmul( orientation, AD%u(1)%HubMotion%Orientation( :,:,1) )
and I wonder if this is the right approach or is there something else to consider?
Could you explain what is the purpose of the following lines of code? I am not fully clear on what is being done here in terms of mapping coordinates or why it is necessary:

BladePtLoads mesh (I am assuming this is in a global coordinate system) is reshuffled into some other coordinate system for what exact reason?
else
RtHSdat%FSAero(1,K,J) = u%BladePtLoads(K)%Force(1,NodeNum) / p%DRNodes(J)
RtHSdat%FSAero(2,K,J) = u%BladePtLoads(K)%Force(3,NodeNum) / p%DRNodes(J)
RtHSdat%FSAero(3,K,J) = u%BladePtLoads(K)%Force(2,NodeNum) / p%DRNodes(J)
RtHSdat%MMAero(1,K,J) = u%BladePtLoads(K)%Moment(1,NodeNum) / p%DRNodes(J)
RtHSdat%MMAero(2,K,J) = u%BladePtLoads(K)%Moment(3,NodeNum) / p%DRNodes(J)
RtHSdat%MMAero(3,K,J) = u%BladePtLoads(K)%Moment(2,NodeNum) / p%DRNodes(J)
end if

Similar swapping is done for the outputs required for AeroDyn but the negative value happen in the second component. Is this mapping not analogous to the one in the CalculateForceMoments routine above?
else
!Translational Displacements
y%BladeLn2Mesh(K)%TranslationDisp(1,NodeNum) = m%RtHS%rS (1,K,J2) ! = the distance from the undeflected tower centerline to the current blade node in the xi ( z1) direction
y%BladeLn2Mesh(K)%TranslationDisp(2,NodeNum) = 1.*m%RtHS%rS (3,K,J2) ! = the distance from the undeflected tower centerline to the current blade node in the yi (z3) direction
y%BladeLn2Mesh(K)%TranslationDisp(3,NodeNum) = m%RtHS%rS (2,K,J2) + p%PtfmRefzt ! = the distance from the nominal tower base position (i.e., the undeflected position of the tower base) to the current blade node in the zi ( z2) direction
!Orientation
y%BladeLn2Mesh(K)%Orientation(1,1,NodeNum) = m%CoordSys%n1(K,J2,1)
y%BladeLn2Mesh(K)%Orientation(2,1,NodeNum) = m%CoordSys%n2(K,J2,1)
y%BladeLn2Mesh(K)%Orientation(3,1,NodeNum) = m%CoordSys%n3(K,J2,1)
y%BladeLn2Mesh(K)%Orientation(1,2,NodeNum) = 1.*m%CoordSys%n1(K,J2,3)
y%BladeLn2Mesh(K)%Orientation(2,2,NodeNum) = 1.*m%CoordSys%n2(K,J2,3)
y%BladeLn2Mesh(K)%Orientation(3,2,NodeNum) = 1.*m%CoordSys%n3(K,J2,3)
y%BladeLn2Mesh(K)%Orientation(1,3,NodeNum) = m%CoordSys%n1(K,J2,2)
y%BladeLn2Mesh(K)%Orientation(2,3,NodeNum) = m%CoordSys%n2(K,J2,2)
y%BladeLn2Mesh(K)%Orientation(3,3,NodeNum) = m%CoordSys%n3(K,J2,2)
end if
Thank you very much for clarifying all these questions.
Best regards,
Juan.