Excerpted from an email correspondence from Serena Thevenin-Illig (LEGOS/IRD):

I just would like to draw your attention on some lines in the code of 
mtm-svd-recon.f where you average some angles in order to estimate the 
averaged phase relatively the reference position (line 845), when doing an 
evolving frequency analysis. One can not average angle because it is not a 
continuous function. (-180+180=0 and should be 180 or -180), we better 
compute the average of the position.

Here are the lines I modified:
I define the averages positions and initialize them to zero.

real *8 xave(mmax,mmax),yave(mmax,mmax)
do i=1,iabv
   do j=1,nmode
         xave(i,j) = 0.d0
         yave(i,j) = 0.d0
   enddo
enddo

Around line 845, The I compute cathesian values of complex U with unit amplitude
           x_pos=cos(atan2(dimag(U(i,j)),dreal(U(i,j))))
           y_pos=sin(atan2(dimag(U(i,j)),dreal(U(i,j))))

I rotate the vector by -1*angle_ref
           x_prim=x_pos*cos(-1*angle_ref)-y_pos*sin(-1*angle_ref)
           y_prim=x_pos*sin(-1*angle_ref)+y_pos*cos(-1*angle_ref)
          
I average the positions, (amplitude remains unity)
           xave(i,j)=xave(i,j)+x_prim
           yave(i,j)=yave(i,j)+y_prim

Around line 1158:
            xave(i,j)=xave(i,j)/float(niter)
            yave(i,j)=yave(i,j)/float(niter)
            avephase(i,j)=180.*atan2(yave(i,j),xave(i,j))/pi

Cordially,
                                                                       .
Serena THEVENIN-ILLIG