diff --git a/Source/WindFarmParametrization/ERF_InitWindFarm.cpp b/Source/WindFarmParametrization/ERF_InitWindFarm.cpp
index 99a6f2f20..6f48a711b 100644
--- a/Source/WindFarmParametrization/ERF_InitWindFarm.cpp
+++ b/Source/WindFarmParametrization/ERF_InitWindFarm.cpp
@@ -125,14 +125,6 @@ WindFarm::init_windfarm_lat_lon (const std::string windfarm_loc_table,
         xloc[it] = xloc[it] - xloc_min + windfarm_x_shift;
         yloc[it] = yloc[it] - yloc_min + windfarm_y_shift;
     }
-
-    FILE* file_xy_loc;
-    file_xy_loc = fopen("file_xy_loc_KingPlains.txt","w");
-
-    for(int it = 0;it<xloc.size(); it++){
-        fprintf(file_xy_loc,"%0.15g %0.15g %0.15g\n", xloc[it], yloc[it], 89.0);
-    }
-    fclose(file_xy_loc);
 }
 
 void
@@ -456,6 +448,8 @@ WindFarm::fill_SMark_multifab(const Geometry& geom,
 
             Real z = ProbLoArr[2] + (kk+0.5) * dx[2];
 
+            int turb_indices_overlap[2];
+            int check_int = 0;
             for(int it=0; it<num_turb; it++){
                 Real x0 = d_xloc_ptr[it] + d_sampling_distance*nx;
                 Real y0 = d_yloc_ptr[it] + d_sampling_distance*ny;
@@ -467,13 +461,21 @@ WindFarm::fill_SMark_multifab(const Geometry& geom,
                 }
                 x0 = d_xloc_ptr[it];
                 y0 = d_yloc_ptr[it];
+                //printf("Values are %d, %0.15g, %0.15g\n", it, x0, y0);
 
                 is_cell_marked = find_if_marked(x1, x2, y1, y2, x0, y0,
                                                 nx, ny, d_hub_height, d_rotor_rad, z);
                 if(is_cell_marked) {
                     SMark_array(i,j,k,1) = it;
+                    turb_indices_overlap[check_int] = it;
+                    check_int++;
+                    if(check_int > 1){
+                        printf("Actuator disks with indices %d and %d are overlapping\n",
+                               turb_indices_overlap[0],turb_indices_overlap[1]);
+                        amrex::Error("Actuator disks are overlapping. Visualize actuator_disks.vtk "
+                        " and check the windturbine locations input file. Exiting..");
+                    }
                 }
-
             }
         });
     }
diff --git a/Source/WindFarmParametrization/GeneralActuatorDisk/ERF_AdvanceGeneralAD.cpp b/Source/WindFarmParametrization/GeneralActuatorDisk/ERF_AdvanceGeneralAD.cpp
index 42f9d4be3..aaa64ceb2 100644
--- a/Source/WindFarmParametrization/GeneralActuatorDisk/ERF_AdvanceGeneralAD.cpp
+++ b/Source/WindFarmParametrization/GeneralActuatorDisk/ERF_AdvanceGeneralAD.cpp
@@ -239,7 +239,7 @@ compute_source_terms_Fn_Ft (const Real rad,
         AMREX_ALWAYS_ASSERT(std::fabs(std::exp(-fhub))<=1.0);
         AMREX_ALWAYS_ASSERT(std::fabs(std::exp(-ftip))<=1.0);
 
-        F = 1.0;//2.0/PI*(std::acos(std::exp(-ftip)) + std::acos(std::exp(-fhub)) );
+        F = 2.0/PI*(std::acos(std::exp(-ftip)) + std::acos(std::exp(-fhub)) );
 
         at_new = 1.0/ ( 4.0*F*std::sin(psi)*std::cos(psi)/(s*Ct+1e-10) - 1.0 );
         an_new = 1.0/ ( 1.0 + 4.0*F*std::pow(std::sin(psi),2)/(s*Cn + 1e-10) );
@@ -417,7 +417,6 @@ GeneralAD::source_terms_cellcentered (const Geometry& geom,
             Real y   = ProbLoArr[1] + (jj+0.5)*dx[1];
             Real z   = ProbLoArr[2] + (kk+0.5)*dx[2];
             // ?? Density needed here
-            Real inv_dens_vol = 1.0/(1.0*dx[0]*dx[1]*dx[2]);
 
             int check_int = 0;
 
@@ -468,17 +467,19 @@ GeneralAD::source_terms_cellcentered (const Geometry& geom,
                                                                d_blade_pitch_ptr,
                                                                n_spec_extra);
 
-                        Real Fn = Fn_and_Ft[0];
-                        Real Ft = Fn_and_Ft[1];
+                        Real Fn = 3.0*Fn_and_Ft[0];
+                        Real Ft = 3.0*Fn_and_Ft[1];
                         // Compute the source terms - pass in radial distance, free stream velocity
 
                         Real Fx = Fn*std::cos(phi) + Ft*std::sin(zeta)*std::sin(phi);
                         Real Fy = Fn*std::sin(phi) - Ft*std::sin(zeta)*std::cos(phi);
                         Real Fz = -Ft*std::cos(zeta);
 
-                        source_x = -Fx*inv_dens_vol;
-                        source_y = -Fy*inv_dens_vol;
-                        source_z = -Fz*inv_dens_vol;
+                        //Real dn = (
+
+                        source_x = -Fx/(2.0*PI*rad*dx[0])*1.0/std::pow(2.0*PI,0.5);
+                        source_y = -Fy/(2.0*PI*rad*dx[0])*1.0/std::pow(2.0*PI,0.5);
+                        source_z = -Fz/(2.0*PI*rad*dx[0])*1.0/std::pow(2.0*PI,0.5);
 
 
                         //printf("Val source_x, is %0.15g, %0.15g, %0.15g %0.15g %0.15g %0.15g\n", rad, Fn, Ft, source_x, source_y, source_z);
diff --git a/Source/WindFarmParametrization/SimpleActuatorDisk/ERF_AdvanceSimpleAD.cpp b/Source/WindFarmParametrization/SimpleActuatorDisk/ERF_AdvanceSimpleAD.cpp
index 38a10d65e..7709f6355 100644
--- a/Source/WindFarmParametrization/SimpleActuatorDisk/ERF_AdvanceSimpleAD.cpp
+++ b/Source/WindFarmParametrization/SimpleActuatorDisk/ERF_AdvanceSimpleAD.cpp
@@ -114,7 +114,7 @@ void SimpleAD::compute_freestream_velocity(const MultiFab& cons_in,
                                  amrex::ParallelContext::CommunicatorAll());
 
     get_turb_loc(xloc, yloc);
-    if (ParallelDescriptor::IOProcessor()){
+    /*if (ParallelDescriptor::IOProcessor()){
         for(int it=0; it<xloc.size(); it++){
 
             std::cout << "turbine index, freestream velocity is " << it << " " << freestream_velocity[it] << " " <<
@@ -122,7 +122,7 @@ void SimpleAD::compute_freestream_velocity(const MultiFab& cons_in,
                                                                 freestream_velocity[it]/(disk_cell_count[it] + 1e-10) << " " <<
                                                                 freestream_phi[it]/(disk_cell_count[it] + 1e-10) << "\n";
         }
-    }
+    }*/
 }
 
 void
@@ -195,12 +195,13 @@ SimpleAD::source_terms_cellcentered (const Geometry& geom,
             int jj = amrex::min(amrex::max(j, domlo_y), domhi_y);
             int kk = amrex::min(amrex::max(k, domlo_z), domhi_z);
 
-            int check_int = 0;
 
             Real source_x = 0.0;
             Real source_y = 0.0;
 
-            for(long unsigned int it=0;it<nturbs;it++) {
+            int it = static_cast<int>(SMark_array(ii,jj,kk,1));
+
+              if(it != -1) {
                 Real avg_vel  = d_freestream_velocity_ptr[it]/(d_disk_cell_count_ptr[it] + 1e-10);
                 Real phi      = d_freestream_phi_ptr[it]/(d_disk_cell_count_ptr[it] + 1e-10);
 
@@ -210,8 +211,6 @@ SimpleAD::source_terms_cellcentered (const Geometry& geom,
                     a = 0.5 - 0.5*std::pow(1.0-C_T,0.5);
                 }
                 Real Uinfty_dot_nhat = avg_vel*(std::cos(phi)*nx + std::sin(phi)*ny);
-                if(SMark_array(ii,jj,kk,1) == static_cast<double>(it)) {
-                    check_int++;
                     if(C_T <= 1) {
                         source_x = -2.0*std::pow(Uinfty_dot_nhat, 2.0)*a*(1.0-a)*dx[1]*dx[2]*std::cos(d_turb_disk_angle)/(dx[0]*dx[1]*dx[2])*std::cos(phi);
                         source_y = -2.0*std::pow(Uinfty_dot_nhat, 2.0)*a*(1.0-a)*dx[1]*dx[2]*std::cos(d_turb_disk_angle)/(dx[0]*dx[1]*dx[2])*std::sin(phi);
@@ -220,12 +219,7 @@ SimpleAD::source_terms_cellcentered (const Geometry& geom,
                         source_x = -0.5*C_T*std::pow(Uinfty_dot_nhat, 2.0)*dx[1]*dx[2]*std::cos(d_turb_disk_angle)/(dx[0]*dx[1]*dx[2])*std::cos(phi);
                         source_y = -0.5*C_T*std::pow(Uinfty_dot_nhat, 2.0)*dx[1]*dx[2]*std::cos(d_turb_disk_angle)/(dx[0]*dx[1]*dx[2])*std::sin(phi);
                     }
-                }
-            }
-            if(check_int > 1){
-                amrex::Error("Actuator disks are overlapping. Visualize actuator_disks.vtk "
-                             "and check the windturbine locations input file. Exiting..");
-            }
+             }
 
             simpleAD_array(i,j,k,0) = source_x;
             simpleAD_array(i,j,k,1) = source_y;