pph inmom set. if statement for mix_mode reogranaized

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)
@@ -99,10 +98,11 @@ module obl_inmom
       !print *, "neutral_mld:", neutral_mld
     end if
     
+    ! 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
-          if (kh_km_mode == 1) 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,:))
+          end if
+        end do
+      end do
+    ! obl_pph mixing mode
     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)
+      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)
+          end if
+        end do
+      end do
+      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 if
         end do
       end do
-
-    if (kh_km_mode == 2) then
-      kh = kh * 10000.0
-      km = km * 10000.0
     end if
     
     ! print *, "Kh first:", kh(3,3,1:4)