diff --git a/obl_inmom.f90 b/obl_inmom.f90
index 83169ab4831f055a0d7e153f43d267a5bee5fe7c..b51725f401bc98af901e20d674d2b249879ef792 100644
--- a/obl_inmom.f90
+++ b/obl_inmom.f90
@@ -76,15 +76,14 @@ module obl_inmom
     kh_b0 = anubgrt
     km_b0 = anubgru
 
+    !! pph_vars
     nx = size(rit, 1)
     ny = size(rit, 2)
     nz = size(rit, 3)
 
-    ! print *, "hello 1"
     den = legacy_denp(tt, ss + 35.0, 0.0)
 
     if (richnum_mode == 1) then
-      ! print *, "hello rit calc"
       call legacy_u2(uu, dy, dyh, hhu, hhq, border_shift, lu, u2)
       call legacy_v2(vv, dx, dxh, hhv, hhq, border_shift, lu, v2)
       call legacy_s2(u2, v2, lu, s2)
@@ -98,11 +97,12 @@ module obl_inmom
       call sync_xy_border_2d(u_dynH)
       !print *, "neutral_mld:", neutral_mld
     end if
-
-    do j = 1, ny
-      do i = 1, nx
-        if (lu(i, j) > lu_min) then
-          if (kh_km_mode == 1) then
+    
+    ! obl_legacy mixing mode
+    if (kh_km_mode == 1) then
+      do j = 1, ny
+        do i = 1, nx
+          if (lu(i, j) > lu_min) then
             call legacy_str("kh", taux(i,j), tauy(i,j), rh0, kh_str)
             call legacy_undimdepth("kh", kh_str, hhq(i,j), aice0(i,j), kh_undimdepth)
             call legacy_kh_b(zw(:), hhq(i,j), kh_unstable, kh_undimdepth, kh_b0, kh_b)
@@ -113,22 +113,37 @@ module obl_inmom
             call legacy_undimdepth("km", km_str, hhq(i,j), aice0(i,j), km_undimdepth)
             call legacy_km_b(zw, km_unstable, km_undimdepth, km_b0, km_b)
             call legacy_km(rit(i,j,:), km_0, km_b(:), km(i,j,:))
-          else if (kh_km_mode == 2) then
-            ! call get_Kh(kh(i,j,:), rit(i,j,:), nz)
-            ! call get_Km(km(i,j,:), rit(i,j,:), nz)
+          end if
+        end do
+      end do
+    ! obl_pph mixing mode
+    else if (kh_km_mode == 2) then
+      pphParams%Km_0 = 7.0 * 0.01
+      pphParams%Kh_0 = 5.0 * 0.01
+      pphParams%alpha = 5.0
+      pphParams%n = 2.0
+      pphParams%Km_b = 5e-6
+      pphParams%Kh_b = 0.00001
+      do j = 1, ny
+        do i = 1, nx
+          if (lu(i, j) > lu_min) then
             call pph_kh(kh(i,j,:), rit(i,j,:), pphParams, nz)
             call pph_km(km(i,j,:), rit(i,j,:), pphParams, nz)
-          else if (kh_km_mode == 3) then
-            call pph_dyn_kh(kh(i,j,:), rit(i,j,:), u_dynH(i,j), neutral_mld(i,j), pphdynParams, nz)
-            call pph_dyn_km(km(i,j,:), rit(i,j,:), u_dynH(i,j), neutral_mld(i,j), pphdynParams, nz)
           end if
-        end if
+        end do
       end do
-    end do
-
-    if (kh_km_mode == 2) then
       kh = kh * 10000.0
       km = km * 10000.0
+    ! obl_pph_dyn mixing mode
+    else if (kh_km_mode == 3) then
+      do j = 1, ny
+        do i = 1, nx
+          if (lu(i, j) > lu_min) then
+            call pph_dyn_kh(kh(i,j,:), rit(i,j,:), u_dynH(i,j), neutral_mld(i,j), pphdynParams, nz)
+            call pph_dyn_km(km(i,j,:), rit(i,j,:), u_dynH(i,j), neutral_mld(i,j), pphdynParams, nz)
+          end if
+        end do
+      end do
     end if
     
     ! print *, "Kh first:", kh(3,3,1:4)