diff --git a/examples/solo_ocean/double_gyre/input.nml b/examples/solo_ocean/double_gyre/input.nml
index b4eced0b15..8d1bd6a9c5 100644
--- a/examples/solo_ocean/double_gyre/input.nml
+++ b/examples/solo_ocean/double_gyre/input.nml
@@ -6,7 +6,4 @@
          parameter_filename = 'MOM_input',
                               'MOM_override' /
 
- &fms_nml
-         domains_stack_size = 710000,
-         stack_size = 0 /
  
diff --git a/src/ALE/MOM_ALE.F90 b/src/ALE/MOM_ALE.F90
index 0e30dfc250..02bb11c457 100644
--- a/src/ALE/MOM_ALE.F90
+++ b/src/ALE/MOM_ALE.F90
@@ -320,6 +320,7 @@ subroutine pressure_gradient_plm( CS, S_t, S_b, T_t, T_b, G, tv, h )
   ! Local variables
   integer :: i, j, k
   real    :: hTmp(G%ke)
+  real    :: tmp(G%ke)
   real, dimension(CS%nk,2) :: &
       ppoly_linear_E            !Edge value of polynomial
   real, dimension(CS%nk,CS%degree_linear+1) :: &
@@ -331,19 +332,18 @@ subroutine pressure_gradient_plm( CS, S_t, S_b, T_t, T_b, G, tv, h )
   ! in 'ALE_memory_allocation'.
 
   ! Determine reconstruction within each column
-!$OMP parallel do default(shared) private(hTmp,ppoly_linear_E,ppoly_linear_coefficients)
-  do i = G%isc,G%iec+1
-    do j = G%jsc,G%jec+1
-     
+!$OMP parallel do default(shared) private(hTmp,ppoly_linear_E,ppoly_linear_coefficients,tmp)
+  do j = G%jsc,G%jec+1
+    do i = G%isc,G%iec+1
       ! Build current grid
       hTmp(:) = h(i,j,:)
-      
+      tmp(:) = tv%S(i,j,:)
       ! Reconstruct salinity profile    
       ppoly_linear_E = 0.0
       ppoly_linear_coefficients = 0.0
-      call PLM_reconstruction( G%ke, hTmp, tv%S(i,j,:), ppoly_linear_E, ppoly_linear_coefficients )
+      call PLM_reconstruction( G%ke, hTmp, tmp, ppoly_linear_E, ppoly_linear_coefficients )
       if (CS%boundary_extrapolation_for_pressure) call &
-        PLM_boundary_extrapolation( G%ke, hTmp, tv%S(i,j,:), ppoly_linear_E, ppoly_linear_coefficients )
+        PLM_boundary_extrapolation( G%ke, hTmp, tmp, ppoly_linear_E, ppoly_linear_coefficients )
       
       do k = 1,G%ke
         S_t(i,j,k) = ppoly_linear_E(k,1)
@@ -353,9 +353,10 @@ subroutine pressure_gradient_plm( CS, S_t, S_b, T_t, T_b, G, tv, h )
       ! Reconstruct temperature profile 
       ppoly_linear_E = 0.0
       ppoly_linear_coefficients = 0.0
-      call PLM_reconstruction( G%ke, hTmp, tv%T(i,j,:), ppoly_linear_E, ppoly_linear_coefficients )
+      tmp(:) = tv%T(i,j,:)
+      call PLM_reconstruction( G%ke, hTmp, tmp, ppoly_linear_E, ppoly_linear_coefficients )
       if (CS%boundary_extrapolation_for_pressure) call &
-        PLM_boundary_extrapolation( G%ke, hTmp, tv%T(i,j,:), ppoly_linear_E, ppoly_linear_coefficients )
+        PLM_boundary_extrapolation( G%ke, hTmp, tmp, ppoly_linear_E, ppoly_linear_coefficients )
       
       do k = 1,G%ke
         T_t(i,j,k) = ppoly_linear_E(k,1)
@@ -393,6 +394,7 @@ subroutine pressure_gradient_ppm( CS, S_t, S_b, T_t, T_b, G, tv, h )
   ! Local variables
   integer :: i, j, k
   real    :: hTmp(G%ke)
+  real    :: tmp(G%ke)
   real, dimension(CS%nk,2) :: &
       ppoly_parab_E            !Edge value of polynomial
   real, dimension(CS%nk,CS%degree_parab+1) :: &
@@ -405,19 +407,20 @@ subroutine pressure_gradient_ppm( CS, S_t, S_b, T_t, T_b, G, tv, h )
 
   ! Determine reconstruction within each column
 !$OMP parallel do default(shared) private(hTmp,ppoly_parab_E,ppoly_parab_coefficients)
-  do i = G%isc,G%iec+1
-    do j = G%jsc,G%jec+1
+  do j = G%jsc,G%jec+1
+    do i = G%isc,G%iec+1
      
       ! Build current grid
       hTmp(:) = h(i,j,:)
+      tmp(:) = tv%S(i,j,:)
       
       ! Reconstruct salinity profile    
       ppoly_parab_E = 0.0
       ppoly_parab_coefficients = 0.0
-      call edge_values_implicit_h4( G%ke, hTmp, tv%S(i,j,:), CS%edgeValueWrk, ppoly_parab_E )
-      call PPM_reconstruction( G%ke, hTmp, tv%S(i,j,:), ppoly_parab_E, ppoly_parab_coefficients )
+      call edge_values_implicit_h4( G%ke, hTmp, tmp, CS%edgeValueWrk, ppoly_parab_E )
+      call PPM_reconstruction( G%ke, hTmp, tmp, ppoly_parab_E, ppoly_parab_coefficients )
       if (CS%boundary_extrapolation_for_pressure) call &
-        PPM_boundary_extrapolation( G%ke, hTmp, tv%S(i,j,:), ppoly_parab_E, ppoly_parab_coefficients )
+        PPM_boundary_extrapolation( G%ke, hTmp, tmp, ppoly_parab_E, ppoly_parab_coefficients )
       
       do k = 1,G%ke
         S_t(i,j,k) = ppoly_parab_E(k,1)
@@ -427,10 +430,11 @@ subroutine pressure_gradient_ppm( CS, S_t, S_b, T_t, T_b, G, tv, h )
       ! Reconstruct temperature profile 
       ppoly_parab_E = 0.0
       ppoly_parab_coefficients = 0.0
-      call edge_values_implicit_h4( G%ke, hTmp, tv%T(i,j,:), CS%edgeValueWrk, ppoly_parab_E )
-      call PPM_reconstruction( G%ke, hTmp, tv%T(i,j,:), ppoly_parab_E, ppoly_parab_coefficients )
+      tmp(:) = tv%T(i,j,:)
+      call edge_values_implicit_h4( G%ke, hTmp, tmp, CS%edgeValueWrk, ppoly_parab_E )
+      call PPM_reconstruction( G%ke, hTmp, tmp, ppoly_parab_E, ppoly_parab_coefficients )
       if (CS%boundary_extrapolation_for_pressure) call &
-        PPM_boundary_extrapolation( G%ke, hTmp, tv%T(i,j,:), ppoly_parab_E, ppoly_parab_coefficients )
+        PPM_boundary_extrapolation( G%ke, hTmp, tmp, ppoly_parab_E, ppoly_parab_coefficients )
       
       do k = 1,G%ke
         T_t(i,j,k) = ppoly_parab_E(k,1)
diff --git a/src/core/MOM.F90 b/src/core/MOM.F90
index 46b3e35cb1..ea3d81133a 100644
--- a/src/core/MOM.F90
+++ b/src/core/MOM.F90
@@ -1257,7 +1257,7 @@ subroutine step_MOM(fluxes, state, Time_start, time_interval, CS)
 
     if (showCallTree) call callTree_leave("DT cycles (step_MOM)")
 
-  call cpu_clock_end(id_clock_other)
+    call cpu_clock_end(id_clock_other)
 
   enddo ! End of n loop
 
@@ -2375,6 +2375,7 @@ subroutine calculate_surface_state(state, u, v, h, ssh, G, CS, p_atm)
 
     depth_ml = CS%Hmix
   !   Determine the mean properties of the uppermost depth_ml fluid.
+!$OMP parallel do default(shared) private(depth,dh)
     do j=js,je
       do i=is,ie
         depth(i) = 0.0
@@ -2410,7 +2411,7 @@ subroutine calculate_surface_state(state, u, v, h, ssh, G, CS, p_atm)
         else
           state%sfc_density(i,j) = state%sfc_density(i,j) / depth(i)
         endif
-        state%Hml(:,:) = depth(i)
+        state%Hml(i,j) = depth(i)
       enddo
     enddo ! end of j loop
   endif                                             ! end BULKMIXEDLAYER
@@ -2421,13 +2422,6 @@ subroutine calculate_surface_state(state, u, v, h, ssh, G, CS, p_atm)
   state%TempxPmE => CS%tv%TempxPmE
   state%internal_heat => CS%tv%internal_heat
 
-  if (associated(state%salt_deficit) .and. associated(CS%tv%salt_deficit)) then
-    do j=js,je ; do i=is,ie
-      ! Convert from gSalt to kgSalt
-      state%salt_deficit(i,j) = 1000.0 * CS%tv%salt_deficit(i,j)
-    enddo ; enddo
-  endif
-
   ! Allocate structures for ocean_mass, ocean_heat, and ocean_salt.  This could
   ! be wrapped in a run-time flag to disable it for economy, since the 3-d
   ! sums are not negligible.
@@ -2443,13 +2437,24 @@ subroutine calculate_surface_state(state, u, v, h, ssh, G, CS, p_atm)
     endif
   endif
 
+!$OMP parallel default(shared) private(mass)
+  if (associated(state%salt_deficit) .and. associated(CS%tv%salt_deficit)) then
+!$OMP do
+    do j=js,je ; do i=is,ie
+      ! Convert from gSalt to kgSalt
+      state%salt_deficit(i,j) = 1000.0 * CS%tv%salt_deficit(i,j)
+    enddo ; enddo
+  endif
+
   if (associated(state%ocean_mass) .and. associated(state%ocean_heat) .and. &
       associated(state%ocean_salt)) then
+!$OMP do
     do j=js,je ; do i=is,ie
       state%ocean_mass(i,j) = 0.0
       state%ocean_heat(i,j) = 0.0 ; state%ocean_salt(i,j) = 0.0
     enddo ; enddo
-    do k=1,nz ; do j=js,je ; do i=is,ie
+!$OMP do
+    do j=js,je ; do k=1,nz; do i=is,ie
       mass = G%H_to_kg_m2*h(i,j,k)
       state%ocean_mass(i,j) = state%ocean_mass(i,j) + mass
       state%ocean_heat(i,j) = state%ocean_heat(i,j) + mass*CS%tv%T(i,j,k)
@@ -2458,27 +2463,34 @@ subroutine calculate_surface_state(state, u, v, h, ssh, G, CS, p_atm)
     enddo ; enddo ; enddo
   else
     if (associated(state%ocean_mass)) then
+!$OMP do
       do j=js,je ; do i=is,ie ; state%ocean_mass(i,j) = 0.0 ; enddo ; enddo
-      do k=1,nz ; do j=js,je ; do i=is,ie
+!$OMP do
+      do j=js,je ; do k=1,nz ; do i=is,ie
         state%ocean_mass(i,j) = state%ocean_mass(i,j) + G%H_to_kg_m2*h(i,j,k)
       enddo ; enddo ; enddo
     endif
     if (associated(state%ocean_heat)) then
+!$OMP do
       do j=js,je ; do i=is,ie ; state%ocean_heat(i,j) = 0.0 ; enddo ; enddo
-      do k=1,nz ; do j=js,je ; do i=is,ie
+!$OMP do
+      do j=js,je ; do k=1,nz ; do i=is,ie
         mass = G%H_to_kg_m2*h(i,j,k)
         state%ocean_heat(i,j) = state%ocean_heat(i,j) + mass*CS%tv%T(i,j,k)
       enddo ; enddo ; enddo
     endif
     if (associated(state%ocean_salt)) then
+!$OMP do
       do j=js,je ; do i=is,ie ; state%ocean_salt(i,j) = 0.0 ; enddo ; enddo
-      do k=1,nz ; do j=js,je ; do i=is,ie
+!$OMP do
+      do j=js,je ; do k=1,nz ; do i=is,ie
         mass = G%H_to_kg_m2*h(i,j,k)
         state%ocean_salt(i,j) = state%ocean_salt(i,j) + &
                                 mass * (1.0e-3*CS%tv%S(i,j,k))
       enddo ; enddo ; enddo
     endif
   endif
+!$OMP end parallel
 
   if (associated(CS%visc%taux_shelf)) state%taux_shelf => CS%visc%taux_shelf
   if (associated(CS%visc%tauy_shelf)) state%tauy_shelf => CS%visc%tauy_shelf
diff --git a/src/core/MOM_PressureForce_Montgomery.F90 b/src/core/MOM_PressureForce_Montgomery.F90
index fca5db4a93..26368acaf6 100644
--- a/src/core/MOM_PressureForce_Montgomery.F90
+++ b/src/core/MOM_PressureForce_Montgomery.F90
@@ -665,6 +665,7 @@ subroutine Set_pbce_Bouss(e, tv, G, g_Earth, Rho0, GFS_scale, pbce, rho_star)
 
   if (use_EOS) then
     if (present(rho_star)) then
+!$OMP parallel do default(shared) private(Ihtot)
       do j=Jsq,Jeq+1
         do i=Isq,Ieq+1
           Ihtot(i) = 1.0 / ((e(i,j,1)-e(i,j,nz+1)) + h_neglect)
@@ -676,6 +677,7 @@ subroutine Set_pbce_Bouss(e, tv, G, g_Earth, Rho0, GFS_scale, pbce, rho_star)
         enddo ; enddo
       enddo ! end of j loop
     else
+!$OMP parallel do default(shared) private(Ihtot,press,rho_in_situ,T_int,S_int,dR_dT,dR_dS)
       do j=Jsq,Jeq+1
         do i=Isq,Ieq+1
           Ihtot(i) = 1.0 / ((e(i,j,1)-e(i,j,nz+1)) + h_neglect)
@@ -704,6 +706,7 @@ subroutine Set_pbce_Bouss(e, tv, G, g_Earth, Rho0, GFS_scale, pbce, rho_star)
       enddo ! end of j loop
     endif
   else ! not use_EOS
+!$OMP parallel do default(shared) private(Ihtot)
     do j=Jsq,Jeq+1
       do i=Isq,Ieq+1
         Ihtot(i) = 1.0 / ((e(i,j,1)-e(i,j,nz+1)) + h_neglect)
diff --git a/src/core/MOM_forcing_type.F90 b/src/core/MOM_forcing_type.F90
index abbfc95e34..2fe8983807 100644
--- a/src/core/MOM_forcing_type.F90
+++ b/src/core/MOM_forcing_type.F90
@@ -461,6 +461,7 @@ subroutine calculateBuoyancyFlux2d(G, fluxes, optics, h, Temp, Salt, tv, buoyanc
 
   netT(:) = 0. ; netS(:) = 0.
 
+!$OMP parallel do default(shared) firstprivate(netT,netS)
   do j = G%jsc, G%jec
     call calculateBuoyancyFlux1d(G, fluxes, optics, h, Temp, Salt, tv, j, buoyancyFlux(:,j,:), netT, netS )
     if (present(netHeatMinusSW)) netHeatMinusSW(:,j) = netT(:)
diff --git a/src/diagnostics/MOM_wave_speed.F90 b/src/diagnostics/MOM_wave_speed.F90
index 4b55692d28..1dcb7ffc8f 100644
--- a/src/diagnostics/MOM_wave_speed.F90
+++ b/src/diagnostics/MOM_wave_speed.F90
@@ -163,7 +163,8 @@ subroutine wave_speed(h, tv, G, cg1, CS, full_halos)
   rescale = 1024.0**4 ; I_rescale = 1.0/rescale
 
   min_h_frac = tol1 / real(nz)
-
+!$OMP parallel do default(private) shared(is,ie,js,je,nz,h,G,min_h_frac,use_EOS,T,S, &
+!$OMP                                     H_to_pres,tv,cg1,g_Rho0,rescale,I_rescale) 
   do j=js,je
     !   First merge very thin layers with the one above (or below if they are
     ! at the top).  This also transposes the row order so that columns can
diff --git a/src/parameterizations/vertical/MOM_KPP.F90 b/src/parameterizations/vertical/MOM_KPP.F90
index a7db4189b3..cbff828036 100644
--- a/src/parameterizations/vertical/MOM_KPP.F90
+++ b/src/parameterizations/vertical/MOM_KPP.F90
@@ -380,6 +380,10 @@ subroutine KPP_calculate(CS, G, h, Temp, Salt, u, v, EOS, uStar, buoyFlux, Kt, K
 
   if (CS%id_Kd_in > 0) call post_data(CS%id_Kd_in, Kt, CS%diag)
 
+!$OMP parallel do default(private) shared(G,CS,EOS,uStar,Temp,Salt,u,v,h,GoRho,buoyFlux, &
+!$OMP                                     const1,nonLocalTransHeat,       &
+!$OMP                                     nonLocalTransScalar,Kt,Ks,Kv) &
+!$OMP                        firstprivate(nonLocalTrans) 
   do j = G%jsc, G%jec
     do i = G%isc, G%iec
       if (G%mask2dT(i,j)==0.) cycle ! Skip calling KPP for land points
@@ -826,7 +830,7 @@ subroutine KPP_applyNonLocalTransport(CS, G, h, nonLocalTrans, surfFlux, dt, sca
   endif
 
   if (diagHeat .or. diagSalt) dSdt(:,:,:) = 0. ! Zero halos for diagnostics ???
-
+!$OMP parallel do default(shared)
   do k = 1, G%ke
     do j = G%jsc, G%jec
       do i = G%isc, G%iec
diff --git a/src/parameterizations/vertical/MOM_diabatic_driver.F90 b/src/parameterizations/vertical/MOM_diabatic_driver.F90
index 284cb54f6b..a4cc1700c7 100644
--- a/src/parameterizations/vertical/MOM_diabatic_driver.F90
+++ b/src/parameterizations/vertical/MOM_diabatic_driver.F90
@@ -164,6 +164,7 @@ module MOM_diabatic_driver
   type(diag_ctrl), pointer :: diag ! A structure that is used to regulate the
                              ! timing of diagnostic output.
   real :: MLDdensityDifference ! Density difference used to determine MLD_user
+  integer :: nsw               ! SW_NBANDS
   integer :: id_dudt_dia = -1, id_dvdt_dia = -1, id_wd = -1
   integer :: id_ea = -1 , id_eb = -1, id_Kd_z = -1, id_Kd_interface = -1
   integer :: id_Tdif_z = -1, id_Tadv_z = -1, id_Sdif_z = -1, id_Sadv_z = -1
@@ -1413,6 +1414,7 @@ subroutine triDiagTS(G, is, ie, js, je, hold, ea, eb, T, S)
   real :: h_tr, b_denom_1
   integer :: i, j, k
 
+!$OMP parallel do default(shared) private(h_tr,b1,d1,c1,b_denom_1)
   do j=js,je
     do i=is,ie
       h_tr = hold(i,j,1) + G%H_subroundoff
@@ -1733,6 +1735,9 @@ subroutine diabatic_driver_init(Time, G, param_file, useALEalgorithm, diag, &
     endif
   endif
 
+  CS%nsw = 0
+  if (ASSOCIATED(CS%optics)) CS%nsw = CS%optics%nbands
+
 end subroutine diabatic_driver_init
 
 subroutine diabatic_driver_end(CS)
@@ -1915,8 +1920,8 @@ subroutine applyBoundaryFluxes(CS, G, dt, fluxes, optics, ea, h, tv)
   real, dimension(SZI_(G)) :: netThickness, netHeat, netSalt, htot, Ttot
   real, dimension(SZI_(G), SZK_(G)) :: h2d, T2d, eps
   integer, dimension(SZI_(G), SZK_(G)) :: ksort
-  real, allocatable, dimension(:,:) :: Pen_SW_bnd
-  real, allocatable, dimension(:,:,:) :: opacityBand
+  real, dimension(max(CS%nsw,1),SZI_(G))         :: Pen_SW_bnd
+  real, dimension(max(CS%nsw,1),SZI_(G),SZK_(G)) :: opacityBand
   logical :: use_riverHeatContent, useCalvingHeatContent
   integer, parameter :: maxGroundings = 5
   integer :: numberOfGroundings, iGround(maxGroundings), jGround(maxGroundings)
@@ -1928,10 +1933,7 @@ subroutine applyBoundaryFluxes(CS, G, dt, fluxes, optics, ea, h, tv)
   ! Skip applying forcing if fluxes%sw is not associated.
   if (.not.ASSOCIATED(fluxes%sw)) return
 
-  nsw = 0
-  if (ASSOCIATED(optics)) nsw = optics%nbands
-  allocate(Pen_SW_bnd(max(nsw,1),SZI_(G)))
-  allocate(opacityBand(max(nsw,1),SZI_(G),SZK_(G)))
+  nsw = CS%nsw
 
   Irho0 = 1.0 / G%Rho0
   I_Cp = 1.0 / fluxes%C_p
@@ -1957,7 +1959,11 @@ subroutine applyBoundaryFluxes(CS, G, dt, fluxes, optics, ea, h, tv)
   if (CS%id_createdH>0) CS%createdH(:,:) = 0.
 
   numberOfGroundings = 0
+!$OMP parallel do default(shared) private(opacityBand,h2d,T2d,eps,htot,netThickness,  &
+!$OMP                                     netHeat,netSalt,Pen_SW_bnd,fractionOfForcing, &
+!$OMP                                     dThickness,dTemp,dSalt,hOld,Ithickness,Ttot)
   do j=js,je ! Work in vertical slices (this is a hold over from the routines called with a j argument)
+    Ttot = 0
     ! Copy state into 2D-slice arrays
     do k=1,nz
       do i=is,ie
@@ -2026,11 +2032,13 @@ subroutine applyBoundaryFluxes(CS, G, dt, fluxes, optics, ea, h, tv)
           endif
         enddo ! k
         if (netThickness(i)/=0.) then ! If anything remains then we grounded out
+!$OMP critical
           numberOfGroundings = numberOfGroundings +1
+!$OMP end critical
           if (numberOfGroundings<=maxGroundings) then
             iGround(numberOfGroundings) = i ! Record i,j location of event for
             jGround(numberOfGroundings) = j ! warning message
-            hGrounding = netThickness(i)
+            hGrounding(numberOfGroundings) = netThickness(i)
           endif
           if (CS%id_createdH>0) CS%createdH(i,j) = CS%createdH(i,j) - netThickness(i)/dt
         endif
@@ -2057,6 +2065,7 @@ subroutine applyBoundaryFluxes(CS, G, dt, fluxes, optics, ea, h, tv)
     enddo
 
   enddo ! j
+
   if (CS%id_createdH > 0) call post_data(CS%id_createdH, CS%createdH, CS%diag)
 
   if (numberOfGroundings>0) then
@@ -2069,9 +2078,6 @@ subroutine applyBoundaryFluxes(CS, G, dt, fluxes, optics, ea, h, tv)
     enddo
   endif
 
-  deallocate(Pen_SW_bnd)
-  deallocate(opacityBand)
-  
 end subroutine applyBoundaryFluxes
 
 end module MOM_diabatic_driver
diff --git a/src/tracer/MOM_tracer_registry.F90 b/src/tracer/MOM_tracer_registry.F90
index 8faa2ad779..95ef6cf074 100644
--- a/src/tracer/MOM_tracer_registry.F90
+++ b/src/tracer/MOM_tracer_registry.F90
@@ -291,16 +291,19 @@ subroutine tracer_vertdiff(h_old, ea, eb, dt, tr, G, &
   is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = G%ke
 
   h_neglect = G%H_subroundoff
-
   sink_dist = 0.0
   if (present(sink_rate)) sink_dist = (dt*sink_rate) * G%m_to_H
+!$OMP parallel default(shared) private(sink,h_minus_dsink,b_denom_1,b1,d1,h_tr,c1)
+!$OMP do
   do j=js,je; do i=is,ie ; sfc_src(i,j) = 0.0 ; btm_src(i,j) = 0.0 ; enddo; enddo
   if (present(sfc_flux)) then
+!$OMP do 
      do j = js, je; do i = is,ie
         sfc_src(i,j) = (sfc_flux(i,j)*dt) * G%kg_m2_to_H
      enddo; enddo
   endif
   if (present(btm_flux)) then
+!$OMP do
      do j = js, je; do i = is,ie
         btm_src(i,j) = (btm_flux(i,j)*dt) * G%kg_m2_to_H
      enddo; enddo
@@ -308,7 +311,7 @@ subroutine tracer_vertdiff(h_old, ea, eb, dt, tr, G, &
 
 
   if (present(sink_rate)) then
-!$OMP parallel do default(shared) private(sink,h_minus_dsink,b_denom_1,b1,d1,h_tr,c1)
+!$OMP do
     do j=js,je
       ! Find the sinking rates at all interfaces, limiting them if necesary
       ! so that the characteristics do not cross within a timestep.
@@ -376,7 +379,7 @@ subroutine tracer_vertdiff(h_old, ea, eb, dt, tr, G, &
       endif ; enddo ; enddo
     enddo
   else
-!$OMP parallel do default(shared) private(h_tr,b_denom_1,b1,d1,c1)
+!$OMP do 
     do j=js,je
       do i=is,ie ; if (G%mask2dT(i,j) > 0.5) then
         h_tr = h_old(i,j,1) + h_neglect
@@ -407,6 +410,8 @@ subroutine tracer_vertdiff(h_old, ea, eb, dt, tr, G, &
     enddo
   endif
 
+!$OMP end parallel
+
 end subroutine tracer_vertdiff
 
 subroutine MOM_tracer_chksum(mesg, Tr, ntr, G)
diff --git a/src/tracer/ideal_age_example.F90 b/src/tracer/ideal_age_example.F90
index 6db64fdca4..aec6fb0bec 100644
--- a/src/tracer/ideal_age_example.F90
+++ b/src/tracer/ideal_age_example.F90
@@ -469,6 +469,7 @@ subroutine ideal_age_tracer_column_physics(h_old, h_new, ea, eb, fluxes, dt, G,
   enddo
   do m=1,CS%ntr ; if (CS%tracer_ages(m) .and. &
                       (year>=CS%tracer_start_year(m))) then
+!$OMP parallel do default(shared)
     do k=CS%nkml+1,nz ; do j=js,je ; do i=is,ie
       CS%tr(i,j,k,m) = CS%tr(i,j,k,m) + G%mask2dT(i,j)*dt*Isecs_per_year
     enddo ; enddo ; enddo