Fixed compiler flags to allow vectorized procedures

Makefile

@@ -170,12 +170,21 @@ fast:
 	  ln -s $(SWIFTEST_HOME)/Makefile.Defines .; \
 	  ln -s $(SWIFTEST_HOME)/Makefile .; \
 	  make fastdir
-
 	cd $(SWIFTEST_HOME)/src/util; \
 	  rm -f Makefile.Defines Makefile; \
 	  ln -s $(SWIFTEST_HOME)/Makefile.Defines .; \
 	  ln -s $(SWIFTEST_HOME)/Makefile .; \
 	  make fastdir
+	cd $(SWIFTEST_HOME)/src/orbel; \
+	  rm -f Makefile.Defines Makefile; \
+	  ln -s $(SWIFTEST_HOME)/Makefile.Defines .; \
+	  ln -s $(SWIFTEST_HOME)/Makefile .; \
+	  make fastdir
+	cd $(SWIFTEST_HOME)/src/drift; \
+	  rm -f Makefile.Defines Makefile; \
+	  ln -s $(SWIFTEST_HOME)/Makefile.Defines .; \
+	  ln -s $(SWIFTEST_HOME)/Makefile .; \
+	  make fastdir
 
 	cd $(SWIFTEST_HOME)/src/orbel; \
 	  rm -f Makefile.Defines Makefile; \
@@ -207,6 +216,12 @@ fast:
 		$(INSTALL_DATA) *.smod $(SWIFTEST_HOME)/include; \
 		rm -f *.o *.smod
 
+	cd $(SWIFTEST_HOME)/src/helio; \
+		$(FORTRAN) $(FFASTFLAGS) -I$(SWIFTEST_HOME)/include -I$(NETCDF_FORTRAN_HOME)/include $(IMKL) -c helio_drift.f90; \
+		$(AR) rv $(SWIFTEST_HOME)/lib/libswiftest.a *.o *.smod; \
+		$(INSTALL_DATA) *.smod $(SWIFTEST_HOME)/include; \
+		rm -f *.o *.smod
+
 libdir:
 	$(FORTRAN) $(FFLAGS) -I$(SWIFTEST_HOME)/include -I$(NETCDF_FORTRAN_HOME)/include $(IMKL) -c *.f90; \
 	$(AR) rv $(SWIFTEST_HOME)/lib/libswiftest.a *.o *.smod; \

Makefile.Defines

@@ -57,7 +57,7 @@ SIMDVEC = -simd -xhost -align all -assume contiguous_assumed_shape -vecabi=cmdta
 PAR = -qopenmp -parallel  
 HEAPARR = -heap-arrays 4194304 
 OPTREPORT = -qopt-report=5 
-IPRODUCTION = -no-wrap-margin -O3 -ipo -qopt-prefetch=0 -sox $(PAR) $(SIMDVEC) #$(HEAPARR)
+IPRODUCTION = -no-wrap-margin -O3 -qopt-prefetch=0 -sox $(PAR) $(SIMDVEC) #$(HEAPARR)
 
 #gfortran flags
 GDEBUG = -g -Og -fbacktrace -fbounds-check -ffree-line-length-none

examples/symba_mars_disk/param.in

@@ -1,6 +1,6 @@
 !Parameter file for the SyMBA-RINGMOONS test
 T0             0.0 
-TSTOP          1200.0
+TSTOP          12000.0
 DT             600.0
 CB_IN          cb.in
 PL_IN          mars.in

src/orbel/orbel.f90

@@ -254,40 +254,24 @@ real(DP) pure function orbel_flon(e,icapn)
       biga =  (-0.5_DP * b + sq)**(1.0_DP / 3.0_DP)
       bigb = -(+0.5_DP * b + sq)**(1.0_DP / 3.0_DP) 
       x = biga + bigb
-      ! write(6,*) 'cubic = ',x**3 +a*x +b
       orbel_flon = x
       ! If capn is VSMALL (or zero) no need to go further than cubic even for
       ! e =1.
-      if( capn < VSMALL) go to 100
-
-      do i = 1,IMAX
-         x2 = x * x
-         f = a0 + x * (a1 + x2 * (a3 + x2 * (a5 + x2 * (a7 + x2 * (a9 + x2 * (a11 + x2))))))
-         fp = b1 + x2 * (b3 + x2 * (b5 + x2 * (b7 + x2 * (b9 + x2 * (b11 + 13 * x2)))))
-         dx = -f / fp
-         !   write(6,*) 'i,dx,x,f : '
-         !   write(6,432) i,dx,x,f
-         432   format(1x,i3,3(2x,1p1e22.15))
-         orbel_flon = x + dx
-         !   if we have converged here there's no point in going on
-            if(abs(dx) <= VSMALL) go to 100
+      if( capn >= VSMALL) then
+         do i = 1,IMAX
+            x2 = x**2
+            f = a0 + x * (a1 + x2 * (a3 + x2 * (a5 + x2 * (a7 + x2 * (a9 + x2 * (a11 + x2))))))
+            fp = b1 + x2 * (b3 + x2 * (b5 + x2 * (b7 + x2 * (b9 + x2 * (b11 + 13 * x2)))))
+            dx = -f / fp
+            orbel_flon = x + dx
+            !   if we have converged here there's no point in going on
+            if(abs(dx) <= VSMALL) exit
             x = orbel_flon
-      end do
-
-      ! abnormal return here - we've gone thru the loop
-      ! imax times without convergence
-      if(iflag == 1) then
-         orbel_flon = -orbel_flon
-         capn = -capn
+         end do
       end if
-      !write(*,*) 'flon : returning without complete convergence'
-      diff = e * sinh(orbel_flon) - orbel_flon - capn
-      !write(*,*) 'n, f, ecc*sinh(f) - f - n : '
-      !write(*,*) capn,orbel_flon,diff
-      return
 
       !  normal return here, but check if capn was originally negative
-      100 if(iflag == 1) then
+      if(iflag == 1) then
          orbel_flon = -orbel_flon
          capn = -capn
       end if