Given B{theta_i} attempt to produce an optimal set of theta values with the same size as the input {theta_i}. The output {theta} set should include all extrema (global and local) and as many other grid points as required to meet the requested number.
The general approach is to: 1. Generate a periodic interpolating spline of B(theta). 2. Find all extrema (essential points). 3. Evaluate the spline on a finer theta mesh, controlled by npadd. 4. Construct a set of B values consisted of all extrema and values on the fine mesh. 5. For each unique B value find all values of theta for which this B appears (setset). 6. Eliminate B values, removing the corresponding values of theta, until we have the required number of theta grid points. We choose which B values to remove based on (user controllable) theta and lambda grid resolution metrics. We strongly bias towards keeping extrema.
Instead of splining and evaluating on a finer mesh we could instead pass a higher resolution B(theta) calculated in geometry when ntheta_geometry > ntheta.
Instead of packing extra theta points in step 3, it might make more sense to finely sample B between the global min and max values?
MRH 15/11/2017 Previous versions of gridgen4_2 produced different theta grids on different machines In this version of gridgen4mod I have corrected this bug by making changes in the following subroutines (in gridgen4_2): get_thetarescorrected this routine so that thetares is 0 when thetares is computed using (theta (i)  theta (j)) and theta(i) equals theta (j) to machine error, and the output should be identically 0. get_lambdares  corrected this routine in the same manner as get_thetares, the output is 0 when it should be identically 0. gg4remove modification of the inequalities used to choose which points to delete from the superset of theta points.
The function of this module is to optimise the choice of theta grid points. It does this by taking the input magnetic field and theta grid (usually generated in the geometry module) and spline fitting the magnetic field to interpolate onto a finemesh theta grid.
By using the metrics defined in get_lambdares and get_thetares the subroutine gg4remove removes points from the finemesh theta grid iteratively until the required number of points remain.The theta grid pair with the lowest value for the metric is deleted at each iteration.The reason for the bug was that these metrics were calculated in a sloppy manner, and the routine for finding the lowest value in the metric did not allow for the case that the metrics were identical up to machine precision:  Because the finemesh theta grid has almost equally spaced points the values for the metrics for the theta grid pairs near theta = 0 were almost exactly equal. Depending on the machine, the metric would be slightly lower for one pair of points than another and therefore the choice of which pair of points to delete 1st from the nonduplicate set was random (actually ~reproducible for the same setup but variable between machines/configs).This behaviour led to different theta grids between different machines.
MRH fixed the above bug by biasing the gg4remove procedure to remove the 1st theta grid pair in array element order (of the sorted theta grid superset array) in the case that the metrics were equal up to machine error. The variable precision_bound controls this behaviour, precision_bound equal 0.0 returns this module to its original behaviour. Finite but small precision_bound is enough to remove machine dependence.
run_name.gridgen.200 contains the debug output from this module To see verbose output ensure that debug_output =.true. and verbose_debug_output =.true.
: This file/module needs tidying/pruning. More consistent usage of tolerance comparisons would help. Testing may be tricky but would be very helpful.
FIXME : Add documentation
Type  Intent  Optional  Attributes  Name  

integer,  intent(in)  ::  n  
integer,  intent(in)  ::  nbmag  
real,  intent(in),  dimension (nbmag)  ::  thetain  
real,  intent(in),  dimension (nbmag)  ::  bmagin  
integer,  intent(in)  ::  npadd  
real,  intent(in)  ::  alknob  
real,  intent(in)  ::  epsknob  
real,  intent(in)  ::  bpknob  
real,  intent(in)  ::  extrknob  
real,  intent(in)  ::  thetamax  
real,  intent(in)  ::  deltaw  
real,  intent(in)  ::  widthw  
real,  intent(in)  ::  tension  
integer,  intent(inout)  ::  ntheta  
integer,  intent(inout)  ::  nlambda  
real,  intent(out),  dimension (ntheta+1)  ::  thetagrid  
real,  intent(out),  dimension (ntheta+1)  ::  bmaggrid  
real,  intent(out),  dimension (nlambda)  ::  alambdagrid 
FIXME : Add documentation
Type  Intent  Optional  Attributes  Name  

integer,  intent(in)  ::  n 
Number of two pi domains included in the input grids 

integer,  intent(in)  ::  nbmag 
Size of the input theta and bmag grids 

real,  intent(in),  dimension (nbmag)  ::  thetain 
Initial theta and bmag grids 

real,  intent(in),  dimension (nbmag)  ::  bmagin 
Initial theta and bmag grids 

integer,  intent(in)  ::  npadd 
Number of additional points to insert between the input theta grid points when determining the starting points 

real,  intent(in)  ::  alknob 
Numerical parameters controlling various choices (see gridgen namelist docs) 

real,  intent(in)  ::  epsknob 
Numerical parameters controlling various choices (see gridgen namelist docs) 

real,  intent(in)  ::  bpknob 
Numerical parameters controlling various choices (see gridgen namelist docs) 

real,  intent(in)  ::  extrknob 
Numerical parameters controlling various choices (see gridgen namelist docs) 

real,  intent(in)  ::  thetamax  
real,  intent(in)  ::  deltaw  
real,  intent(in)  ::  widthw  
real,  intent(in)  ::  tension  
integer,  intent(inout)  ::  ntheta 
On entry ntheta = nbmag  1, on exit ntheta is the number of theta grid points in the final grid. We expect this to be unchanged. On entry nbset is the upper limit on the number of unique B values, i.e. nbset = ntheta / 2 + 1. On exit nbset is the number of unique B values in the output grid. 

integer,  intent(inout)  ::  nbset 
On entry ntheta = nbmag  1, on exit ntheta is the number of theta grid points in the final grid. We expect this to be unchanged. On entry nbset is the upper limit on the number of unique B values, i.e. nbset = ntheta / 2 + 1. On exit nbset is the number of unique B values in the output grid. 

real,  intent(out),  dimension (ntheta+1)  ::  thetagrid 
The optimisied theta grid and corresponding B values 

real,  intent(out),  dimension (ntheta+1)  ::  bmaggrid 
The optimisied theta grid and corresponding B values 

real,  intent(out),  dimension (nbset)  ::  bset 
The ordered set of unique B values on the optimised grid 