Overwrite axis with average axis by center of closed line of the boundary in each poloidal plane
SUBROUTINE InitAverageAxis() ! MODULES USE MODgvec_MHD3D_Vars , ONLY:X1_base,X1_a,X1_b USE MODgvec_MHD3D_Vars , ONLY:X2_base,X2_a,X2_b USE MODgvec_MHD3D_Vars , ONLY:average_axis_move IMPLICIT NONE !----------------------------------------------------------------------------------------------------------------------------------- ! INPUT Variables !----------------------------------------------------------------------------------------------------------------------------------- ! LOCAL Variables REAL(wp) :: X1_b_IP(X1_base%f%mn_nyq(1),X1_base%f%mn_nyq(2)) REAL(wp) :: X2_b_IP(X2_base%f%mn_nyq(1),X2_base%f%mn_nyq(2)) INTEGER :: i_m,i_n REAL(wp) :: dl,lint,x1int,x2int !----------------------------------------------------------------------------------------------------------------------------------- ASSOCIATE(m_nyq=>X1_base%f%mn_nyq(1),n_nyq=>X1_base%f%mn_nyq(2)) X1_b_IP(:,:) = RESHAPE(X1_base%f%evalDOF_IP(0,X1_b),(/m_nyq,n_nyq/)) X2_b_IP(:,:) = RESHAPE(X2_base%f%evalDOF_IP(0,X2_b),(/m_nyq,n_nyq/)) DO i_n=1,n_nyq dl=X1_b_IP(m_nyq,i_n)*X2_b_IP(1,i_n)-X1_b_IP(1,i_n)*X2_b_IP(m_nyq,i_n) lint=dl x1int=(X1_b_IP(m_nyq,i_n)+X1_b_IP(1,i_n))*dl x2int=(X2_b_IP(m_nyq,i_n)+X2_b_IP(1,i_n))*dl DO i_m=2,m_nyq dl=SQRT((X1_b_IP(i_m,i_n)-X1_b_IP(i_m-1,i_n))**2+(X2_b_IP(i_m,i_n)-X2_b_IP(i_m-1,i_n))**2) dl=X1_b_IP(i_m-1,i_n)*X2_b_IP(i_m,i_n)-X1_b_IP(i_m,i_n)*X2_b_IP(i_m-1,i_n) lint=lint+dl x1int=x1int+(X1_b_IP(i_m-1,i_n)+X1_b_IP(i_m,i_n))*dl x2int=x2int+(X2_b_IP(i_m-1,i_n)+X2_b_IP(i_m,i_n))*dl END DO ! c_x= 1/(6A) sum_i (x_i-1+x_i)*(x_i-1*y_i - x_i*y_i-1), A=1/2 sum_i (x_i-1*y_i - x_i*y_i-1) X1_b_IP(:,i_n) = x1int/(3.0_wp*lint) + average_axis_move(1) X2_b_IP(:,i_n) = x2int/(3.0_wp*lint) + average_axis_move(2) END DO X1_a = X1_base%f%initDOF(RESHAPE(X1_b_IP,(/X1_base%f%mn_IP/))) X2_a = X2_base%f%initDOF(RESHAPE(X2_b_IP,(/X2_base%f%mn_IP/))) END ASSOCIATE END SUBROUTINE InitAverageAxis