Add ANN-based mesoscale streamfunction parameterization

src/core/MOM.F90

@@ -1233,7 +1233,7 @@ subroutine step_MOM_dynamics(forces, p_surf_begin, p_surf_end, dt, dt_tr_adv, &
         call calc_slope_functions(h, CS%tv, dt, G, GV, US, CS%VarMix, OBC=CS%OBC)
       call thickness_diffuse(h, CS%uhtr, CS%vhtr, CS%tv, dt_tr_adv, G, GV, US, &
                              CS%MEKE, CS%VarMix, CS%CDp, CS%thickness_diffuse_CSp, &
-                             CS%stoch_CS)
+                             CS%stoch_CS, u, v)
       call cpu_clock_end(id_clock_thick_diff)
       call pass_var(h, G%Domain, clock=id_clock_pass, halo=CS%dyn_h_stencil)
       if (showCallTree) call callTree_waypoint("finished thickness_diffuse_first (step_MOM)")
@@ -1384,8 +1384,8 @@ subroutine step_MOM_dynamics(forces, p_surf_begin, p_surf_end, dt, dt_tr_adv, &
       if (CS%VarMix%use_variable_mixing) &
         call calc_slope_functions(h, CS%tv, dt, G, GV, US, CS%VarMix, OBC=CS%OBC)
       call thickness_diffuse(h, CS%uhtr, CS%vhtr, CS%tv, dt, G, GV, US, &
-                             CS%MEKE, CS%VarMix, CS%CDp, CS%thickness_diffuse_CSp, CS%stoch_CS)
-
+                             CS%MEKE, CS%VarMix, CS%CDp, CS%thickness_diffuse_CSp, &
+                             CS%stoch_CS, u, v)
       call cpu_clock_end(id_clock_thick_diff)
       call pass_var(h, G%Domain, clock=id_clock_pass, halo=CS%dyn_h_stencil)
       if (CS%debug) call hchksum(h,"Post-thickness_diffuse h", G%HI, haloshift=1, unscale=GV%H_to_MKS)

src/core/MOM_isopycnal_slopes.F90

@@ -31,7 +31,8 @@ module MOM_isopycnal_slopes
 !> Calculate isopycnal slopes, and optionally return other stratification dependent functions such as N^2
 !! and dz*S^2*g-prime used, or calculable from factors used, during the calculation.
 subroutine calc_isoneutral_slopes(G, GV, US, h, e, tv, dt_kappa_smooth, use_stanley, slope_x, slope_y, &
-                                  N2_u, N2_v, dzu, dzv, dzSxN, dzSyN, halo, OBC, OBC_N2)
+                                  N2_u, N2_v, dzu, dzv, dzSxN, dzSyN, halo, OBC, OBC_N2, &
+                                  drdx_u, drdy_v, drdz_u, drdz_v)
   type(ocean_grid_type),                       intent(in)    :: G    !< The ocean's grid structure
   type(verticalGrid_type),                     intent(in)    :: GV   !< The ocean's vertical grid structure
   type(unit_scale_type),                       intent(in)    :: US   !< A dimensional unit scaling type
@@ -66,6 +67,22 @@ subroutine calc_isoneutral_slopes(G, GV, US, h, e, tv, dt_kappa_smooth, use_stan
   logical,                           optional, intent(in)    :: OBC_N2 !< If present and true, use interior data
                                                                      !! to calculate stratification at open boundary
                                                                      !! condition faces.
+  real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)+1), &
+                                     optional, intent(inout) :: drdx_u !< Zonal density gradient at u
+                                                                       !! along surfaces of constant z
+                                                                       !! (not along isopycnals or
+                                                                       !! model interfaces) [R L-1 ~> kg m-4]
+  real, dimension(SZI_(G),SZJB_(G),SZK_(GV)+1), &
+                                     optional, intent(inout) :: drdy_v !< Meridional density gradient at v
+                                                                       !! along surfaces of constant z
+                                                                       !! (not along isopycnals or
+                                                                       !! model interfaces) [R L-1 ~> kg m-4]
+  real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)+1), &
+                                     optional, intent(inout) :: drdz_u !< Vertical density gradient
+                                                                       !! at u [R Z-1 ~> kg m-4]
+  real, dimension(SZI_(G),SZJB_(G),SZK_(GV)+1), &
+                                     optional, intent(inout) :: drdz_v !< Vertical density gradient
+                                                                       !! at v [R Z-1 ~> kg m-4]
 
   ! Local variables
   real, dimension(SZI_(G), SZJ_(G), SZK_(GV)) :: &
@@ -131,6 +148,7 @@ subroutine calc_isoneutral_slopes(G, GV, US, h, e, tv, dt_kappa_smooth, use_stan
   real :: G_Rho0        ! The gravitational acceleration divided by density [L2 Z-1 T-2 R-1 ~> m4 s-2 kg-1]
 
   logical :: present_N2_u, present_N2_v
+  logical :: present_drdx_u, present_drdy_v
   logical :: local_open_u_BC, local_open_v_BC ! True if u- or v-face OBCs exist anywhere in the global domain.
   logical :: OBC_friendly  ! If true, open boundary conditions are in use and only interior data should
                         ! be used to calculate N2 at OBC faces.
@@ -170,6 +188,8 @@ subroutine calc_isoneutral_slopes(G, GV, US, h, e, tv, dt_kappa_smooth, use_stan
 
   use_EOS = associated(tv%eqn_of_state)
 
+  present_drdx_u = PRESENT(drdx_u)
+  present_drdy_v = PRESENT(drdy_v)
   present_N2_u = PRESENT(N2_u)
   present_N2_v = PRESENT(N2_v)
   G_Rho0 = GV%g_Earth / GV%Rho0
@@ -209,6 +229,24 @@ subroutine calc_isoneutral_slopes(G, GV, US, h, e, tv, dt_kappa_smooth, use_stan
       dzSyN(i,J,nz+1) = 0.
     enddo ; enddo
   endif
+  ! Set boundary values to zero, since they will be at places
+  ! where streamfunction would be zero.
+  if (present_drdx_u) then
+    do j=js,je ; do I=is-1,ie
+      drdx_u(I,j,1) = 0.
+      drdx_u(I,j,nz+1) = 0.
+      drdz_u(I,j,1) = 0.
+      drdz_u(I,j,nz+1) = 0.
+    enddo ; enddo
+  endif
+  if (present_drdy_v) then
+    do J=js-1,je ; do i=is,ie
+      drdy_v(i,J,1) = 0.
+      drdy_v(i,J,nz+1) = 0.
+      drdz_v(i,J,1) = 0.
+      drdz_v(i,J,nz+1) = 0.
+    enddo ; enddo
+  endif
 
   if (use_EOS) then
     if (present(halo)) then
@@ -254,6 +292,7 @@ subroutine calc_isoneutral_slopes(G, GV, US, h, e, tv, dt_kappa_smooth, use_stan
   !$OMP parallel do default(none) shared(nz,is,ie,js,je,IsdB,use_EOS,G,GV,US,pres,T,S,tv,h,e, &
   !$OMP                                  h_neglect,dz_neglect,h_neglect2, &
   !$OMP                                  present_N2_u,G_Rho0,N2_u,slope_x,dzSxN,EOSdom_u,EOSdom_h1, &
+  !$OMP                                  present_drdx_u,drdx_u,drdz_u, &
   !$OMP                                  local_open_u_BC,dzu,OBC,use_stanley,OBC_friendly) &
   !$OMP                          private(drdiA,drdiB,drdkL,drdkR,pres_u,T_u,S_u,      &
   !$OMP                                  drho_dT_u,drho_dS_u,hg2A,hg2B,hg2L,hg2R,haA, &
@@ -416,6 +455,10 @@ subroutine calc_isoneutral_slopes(G, GV, US, h, e, tv, dt_kappa_smooth, use_stan
       endif
 
       slope_x(I,j,K) = slope
+      if (present_drdx_u) then
+        drdx_u(I,j,K) = drdx
+        drdz_u(I,j,K) = drdz
+      endif
       if (present(dzSxN)) &
         dzSxN(I,j,K) = sqrt( GxSpV_u(I) * max(0., (wtL * ( dzaL * drdkL )) &
                                                 + (wtR * ( dzaR * drdkR ))) / (wtL + wtR) ) & ! dz * N
@@ -431,6 +474,7 @@ subroutine calc_isoneutral_slopes(G, GV, US, h, e, tv, dt_kappa_smooth, use_stan
   !$OMP parallel do default(none) shared(nz,is,ie,js,je,IsdB,use_EOS,G,GV,US,pres,T,S,tv, &
   !$OMP                                  h,h_neglect,e,dz_neglect, &
   !$OMP                                  h_neglect2,present_N2_v,G_Rho0,N2_v,slope_y,dzSyN,EOSdom_v, &
+  !$OMP                                  present_drdy_v,drdy_v,drdz_v, &
   !$OMP                                  dzv,local_open_v_BC,OBC,use_stanley,OBC_friendly) &
   !$OMP                          private(drdjA,drdjB,drdkL,drdkR,pres_v,T_v,S_v,      &
   !$OMP                                  drho_dT_v,drho_dS_v,hg2A,hg2B,hg2L,hg2R,haA, &
@@ -599,6 +643,10 @@ subroutine calc_isoneutral_slopes(G, GV, US, h, e, tv, dt_kappa_smooth, use_stan
         slope = slope * max(G%mask2dT(i,j), G%mask2dT(i,j+1))
       endif
       slope_y(i,J,K) = slope
+      if (present_drdy_v) then
+        drdy_v(i,J,K) = drdy
+        drdz_v(i,J,K) = drdz
+      endif
       if (present(dzSyN)) &
         dzSyN(i,J,K) = sqrt( GxSpV_v(i) * max(0., (wtL * ( dzaL * drdkL )) &
                                                 + (wtR * ( dzaR * drdkR ))) / (wtL + wtR) ) & ! dz * N