Complex Fortran Simulation/Estimation Project (Economics)

I have a project that, while it is not physics per se, could use somebody with that particular skill set when it comes to scientific computing. Of course this is not restricted to physicists. I am running into major robustness issues with my code that probably just needs a fresh set of eyes. The final output of this project would be a successful run of the program, although I imagine there could be a series of checkpoint payments along the way. I have attached the code involved. This is an economics project that estimates the parameters of a model of the labor market. The basic idea is as follows: 1. N workers start at t=0 with a random variable drawn from a parameterized distribution which is treated as the individual's ability. 2. Given this distribution of worker ability and model parameters corresponding to how this ability translates into output, firms which differ in their productivity (and therefore how much they value a given worker) open up vacancies at a cost until the expected profit of a new job opening is zero. 3. Workers are randomly selected with some probability (another set of parameters) and matched with jobs according to the distribution determined in step 2. 4. Each time period, workers transition to a new state. Workers can get a raise, change jobs, get fired for being revealed to be of low ability, or otherwise find themselves unemployed. With some probability they "die" and are replaced with fresh workers. There are state variables that dictate these transition probabilities, but an easy solution to the model has eluded me due to the complication that the firms add to the mix. 5. Given these transitions, worker states are updated and firms once again (each time period) open job vacancies given the distribution of worker states to the point that all firms have a zero expected value on an additional created vacancy. These decisions are not independent however and I have had to use gradient based methods to solve for the distribution of jobs. This has been the source of most of my problems. 6. This process iterates until convergence of the distribution of worker states and thus the distribution of firm vacancies. At this point I take the distribution of firms and begin a second simulation. 7. I simulate the career for N workers given the steady state distribution of jobs obtained from steps 1-6. I look at aggregated variables like the wage distribution at varying points in time, average unemployment durations, etc. These variables are weighted and compared to similar characteristics in the real data. 8. A loss function is used to quantify the distance between the simulated data and the actual data. 9. Original parameter values are updated and the entire process from step 1 is restarted. The process ends when the loss is sufficiently small. _________________________________________________________________ What I need: This code is currently running but the optimization is sloppy (I'm an economist) and needs to be more robust to a greater range of parameter values. I need to be able to get a set of parameter values minimizing the simplified loss function I have specified right now. Right now I am just getting errors for many parameter combinations because the algorithm I wrote to solve for the job vacancy distribution is simply not working. You may continue with what I have or if there is something obvious I have missed and the steady state distribution can be solved for in a simpler way I would be happy to hear about that as well. I have attached all code for the simulation as well as an earlier paper that exposits the model if it would help your understanding of the model. The primary code is structured as follows: [login to view URL] - initial parameters and bounds, calls solution algorithm [login to view URL] - parallel accelerated random search algorithm over parameter values minimizing the loss function [login to view URL] - contains all modules relating to the model solution, simulations, and loss function