Merge pull request #1329 from herrwang0/fix-wavedrag-initRayleigh

Fix the bug that Rayleigh_[uv] was not initialized in MOM_barotropic
mom-ocean · Feb 25, 2021 · f06669e · f06669e
2 parents d8733fe + ac88645
commit f06669e
Showing 1 changed file with 62 additions and 20 deletions.
diff --git a/src/core/MOM_barotropic.F90 b/src/core/MOM_barotropic.F90
@@ -955,6 +955,17 @@ subroutine btstep(U_in, V_in, eta_in, dt, bc_accel_u, bc_accel_v, forces, pbce,
     Datv(i,J) = 0.0 ; bt_rem_v(i,J) = 0.0 ; vhbt0(i,J) = 0.0
   enddo ; enddo
 
+  if (CS%linear_wave_drag) then
+    !$OMP parallel do default(shared)
+    do j=CS%jsdw,CS%jedw ; do I=CS%isdw-1,CS%iedw
+      Rayleigh_u(I,j) = 0.0
+    enddo ; enddo
+    !$OMP parallel do default(shared)
+    do J=CS%jsdw-1,CS%jedw ; do i=CS%isdw,CS%iedw
+      Rayleigh_v(i,J) = 0.0
+    enddo ; enddo
+  endif
+
   ! Copy input arrays into their wide-halo counterparts.
   if (interp_eta_PF) then
     !$OMP parallel do default(shared)
@@ -1900,20 +1911,27 @@ subroutine btstep(U_in, V_in, eta_in, dt, bc_accel_u, bc_accel_v, forces, pbce,
         endif ; enddo ; enddo
         !$OMP end do nowait
       endif
+
       !$OMP do schedule(static)
       do J=jsv-1,jev ; do i=isv-1,iev+1
         vel_prev = vbt(i,J)
         vbt(i,J) = bt_rem_v(i,J) * (vbt(i,J) + &
                     dtbt * ((BT_force_v(i,J) + Cor_v(i,J)) + PFv(i,J)))
         vbt_trans(i,J) = trans_wt1*vbt(i,J) + trans_wt2*vel_prev
+      enddo ; enddo
 
-        if (CS%linear_wave_drag) then
+      if (CS%linear_wave_drag) then
+        !$OMP do schedule(static)
+        do J=jsv-1,jev ; do i=isv-1,iev+1
           v_accel_bt(i,J) = v_accel_bt(i,J) + wt_accel(n) * &
               ((Cor_v(i,J) + PFv(i,J)) - vbt(i,J)*Rayleigh_v(i,J))
-        else
+        enddo ; enddo
+      else
+        !$OMP do schedule(static)
+        do J=jsv-1,jev ; do i=isv-1,iev+1
           v_accel_bt(i,J) = v_accel_bt(i,J) + wt_accel(n) * (Cor_v(i,J) + PFv(i,J))
-        endif
-      enddo ; enddo
+        enddo ; enddo
+      endif
 
       if (integral_BT_cont) then
         !$OMP do schedule(static)
@@ -1969,22 +1987,31 @@ subroutine btstep(U_in, V_in, eta_in, dt, bc_accel_u, bc_accel_v, forces, pbce,
         endif ; enddo ; enddo
         !$OMP end do nowait
       endif
+
       !$OMP do schedule(static)
       do j=jsv,jev ; do I=isv-1,iev
         vel_prev = ubt(I,j)
         ubt(I,j) = bt_rem_u(I,j) * (ubt(I,j) + &
              dtbt * ((BT_force_u(I,j) + Cor_u(I,j)) + PFu(I,j)))
         if (abs(ubt(I,j)) < CS%vel_underflow) ubt(I,j) = 0.0
         ubt_trans(I,j) = trans_wt1*ubt(I,j) + trans_wt2*vel_prev
+      enddo ; enddo
+      !$OMP end do nowait
 
-        if (CS%linear_wave_drag) then
+      if (CS%linear_wave_drag) then
+        !$OMP do schedule(static)
+        do j=jsv,jev ; do I=isv-1,iev
           u_accel_bt(I,j) = u_accel_bt(I,j) + wt_accel(n) * &
-              ((Cor_u(I,j) + PFu(I,j)) - ubt(I,j)*Rayleigh_u(I,j))
-        else
+             ((Cor_u(I,j) + PFu(I,j)) - ubt(I,j)*Rayleigh_u(I,j))
+        enddo ; enddo
+        !$OMP end do nowait
+      else
+        !$OMP do schedule(static)
+        do j=jsv,jev ; do I=isv-1,iev
           u_accel_bt(I,j) = u_accel_bt(I,j) + wt_accel(n) * (Cor_u(I,j) + PFu(I,j))
-        endif
-      enddo ; enddo
-      !$OMP end do nowait
+        enddo ; enddo
+        !$OMP end do nowait
+      endif
 
       if (integral_BT_cont) then
         !$OMP do schedule(static)
@@ -2046,14 +2073,20 @@ subroutine btstep(U_in, V_in, eta_in, dt, bc_accel_u, bc_accel_v, forces, pbce,
              dtbt * ((BT_force_u(I,j) + Cor_u(I,j)) + PFu(I,j)))
         if (abs(ubt(I,j)) < CS%vel_underflow) ubt(I,j) = 0.0
         ubt_trans(I,j) = trans_wt1*ubt(I,j) + trans_wt2*vel_prev
+      enddo ; enddo
 
-        if (CS%linear_wave_drag) then
+      if (CS%linear_wave_drag) then
+        !$OMP do schedule(static)
+        do j=jsv-1,jev+1 ; do I=isv-1,iev
           u_accel_bt(I,j) = u_accel_bt(I,j) + wt_accel(n) * &
               ((Cor_u(I,j) + PFu(I,j)) - ubt(I,j)*Rayleigh_u(I,j))
-        else
+        enddo ; enddo
+      else
+        !$OMP do schedule(static)
+        do j=jsv-1,jev+1 ; do I=isv-1,iev
           u_accel_bt(I,j) = u_accel_bt(I,j) + wt_accel(n) * (Cor_u(I,j) + PFu(I,j))
-        endif
-      enddo ; enddo
+        enddo ; enddo
+      endif
 
       if (integral_BT_cont) then
         !$OMP do schedule(static)
@@ -2128,15 +2161,24 @@ subroutine btstep(U_in, V_in, eta_in, dt, bc_accel_u, bc_accel_v, forces, pbce,
              dtbt * ((BT_force_v(i,J) + Cor_v(i,J)) + PFv(i,J)))
         if (abs(vbt(i,J)) < CS%vel_underflow) vbt(i,J) = 0.0
         vbt_trans(i,J) = trans_wt1*vbt(i,J) + trans_wt2*vel_prev
+      enddo ; enddo
+      !$OMP end do nowait
 
-        if (CS%linear_wave_drag) then
+      if (CS%linear_wave_drag) then
+        !$OMP do schedule(static)
+        do J=jsv-1,jev ; do i=isv,iev
           v_accel_bt(i,J) = v_accel_bt(i,J) + wt_accel(n) * &
-              ((Cor_v(i,J) + PFv(i,J)) - vbt(i,J)*Rayleigh_v(i,J))
-        else
+             ((Cor_v(i,J) + PFv(i,J)) - vbt(i,J)*Rayleigh_v(i,J))
+        enddo ; enddo
+        !$OMP end do nowait
+      else
+        !$OMP do schedule(static)
+        do J=jsv-1,jev ; do i=isv,iev
           v_accel_bt(i,J) = v_accel_bt(i,J) + wt_accel(n) * (Cor_v(i,J) + PFv(i,J))
-        endif
-      enddo ; enddo
-      !$OMP end do nowait
+        enddo ; enddo
+        !$OMP end do nowait
+      endif
+
       if (integral_BT_cont) then
         !$OMP do schedule(static)
         do J=jsv-1,jev ; do i=isv,iev