CPT updates for diagnostics and configurations

src/core/MOM_dynamics_split_RK2.F90

@@ -1219,6 +1219,9 @@ subroutine initialize_dyn_split_RK2(u, v, h, uh, vh, eta, Time, G, GV, US, param
   Accel_diag%PFv   => CS%PFv
   Accel_diag%CAu   => CS%CAu
   Accel_diag%CAv   => CS%CAv
+  Accel_diag%u_accel_bt => CS%u_accel_bt
+  Accel_diag%v_accel_bt => CS%v_accel_bt
+
 
 !  Accel_diag%pbce => CS%pbce
 !  Accel_diag%u_accel_bt => CS%u_accel_bt ; Accel_diag%v_accel_bt => CS%v_accel_bt

src/core/MOM_variables.F90

@@ -171,7 +171,9 @@ module MOM_variables
     du_dt_visc => NULL(), &!< Zonal acceleration due to vertical viscosity [L T-2 ~> m s-2]
     dv_dt_visc => NULL(), &!< Meridional acceleration due to vertical viscosity [L T-2 ~> m s-2]
     du_dt_dia => NULL(), & !< Zonal acceleration due to diapycnal  mixing [L T-2 ~> m s-2]
-    dv_dt_dia => NULL()    !< Meridional acceleration due to diapycnal  mixing [L T-2 ~> m s-2]
+    dv_dt_dia => NULL(), & !< Meridional acceleration due to diapycnal  mixing [L T-2 ~> m s-2]
+    u_accel_bt => NULL(), &!< Pointer to the zonal barotropic-solver acceleration [L T-2 ~> m s-2]
+    v_accel_bt => NULL()   !< Pointer to the meridional barotropic-solver acceleration [L T-2 ~> m s-2]
   real, pointer, dimension(:,:,:) :: du_other => NULL()
                            !< Zonal velocity changes due to any other processes that are
                            !! not due to any explicit accelerations [L T-1 ~> m s-1].

src/diagnostics/MOM_diagnostics.F90

@@ -99,6 +99,7 @@ module MOM_diagnostics
     KE        => NULL(), &  !< KE per unit mass [L2 T-2 ~> m2 s-2]
     dKE_dt    => NULL(), &  !< time derivative of the layer KE [H L2 T-3 ~> m3 s-3]
     PE_to_KE  => NULL(), &  !< potential energy to KE term [m3 s-3]
+    KE_BT     => NULL(), &  !< barotropic contribution to KE term [m3 s-3]
     KE_CorAdv => NULL(), &  !< KE source from the combined Coriolis and
                             !! advection terms [H L2 T-3 ~> m3 s-3].
                             !! The Coriolis source should be zero, but is not due to truncation
@@ -117,7 +118,8 @@ module MOM_diagnostics
   integer :: id_hf_du_dt_2d    = -1, id_hf_dv_dt_2d    = -1
   integer :: id_col_ht         = -1, id_dh_dt          = -1
   integer :: id_KE             = -1, id_dKEdt          = -1
-  integer :: id_PE_to_KE       = -1, id_KE_Coradv      = -1
+  integer :: id_PE_to_KE       = -1, id_KE_BT          = -1
+  integer :: id_KE_Coradv      = -1
   integer :: id_KE_adv         = -1, id_KE_visc        = -1
   integer :: id_KE_horvisc     = -1, id_KE_dia         = -1
   integer :: id_uh_Rlay        = -1, id_vh_Rlay        = -1
@@ -994,9 +996,9 @@ subroutine calculate_energy_diagnostics(u, v, h, uh, vh, ADp, CDp, G, GV, US, CS
   endif
 
   if (.not.G%symmetric) then
-    if (associated(CS%dKE_dt) .OR. associated(CS%PE_to_KE) .OR. associated(CS%KE_CorAdv) .OR. &
-       associated(CS%KE_adv) .OR. associated(CS%KE_visc)  .OR. associated(CS%KE_horvisc).OR. &
-       associated(CS%KE_dia) ) then
+    if (associated(CS%dKE_dt) .OR. associated(CS%PE_to_KE) .OR. associated(CS%KE_BT) .OR. &
+       associated(CS%KE_CorAdv) .OR. associated(CS%KE_adv) .OR. associated(CS%KE_visc) .OR. &
+       associated(CS%KE_horvisc) .OR. associated(CS%KE_dia) ) then
         call create_group_pass(CS%pass_KE_uv, KE_u, KE_v, G%Domain, To_North+To_East)
     endif
   endif
@@ -1040,6 +1042,24 @@ subroutine calculate_energy_diagnostics(u, v, h, uh, vh, ADp, CDp, G, GV, US, CS
     if (CS%id_PE_to_KE > 0) call post_data(CS%id_PE_to_KE, CS%PE_to_KE, CS%diag)
   endif
 
+  if (associated(CS%KE_BT)) then
+    do k=1,nz
+      do j=js,je ; do I=Isq,Ieq
+        KE_u(I,j) = uh(I,j,k) * G%dxCu(I,j) * ADp%u_accel_bt(I,j,k)
+      enddo ; enddo
+      do J=Jsq,Jeq ; do i=is,ie
+        KE_v(i,J) = vh(i,J,k) * G%dyCv(i,J) * ADp%v_accel_bt(i,J,k)
+      enddo ; enddo
+      if (.not.G%symmetric) &
+         call do_group_pass(CS%pass_KE_uv, G%domain)
+      do j=js,je ; do i=is,ie
+        CS%KE_BT(i,j,k) = 0.5 * G%IareaT(i,j) &
+            * (KE_u(I,j) + KE_u(I-1,j) + KE_v(i,J) + KE_v(i,J-1))
+      enddo ; enddo
+    enddo
+    if (CS%id_KE_BT > 0) call post_data(CS%id_KE_BT, CS%KE_BT, CS%diag)
+  endif
+
   if (associated(CS%KE_CorAdv)) then
     do k=1,nz
       do j=js,je ; do I=Isq,Ieq
@@ -1519,6 +1539,7 @@ subroutine MOM_diagnostics_init(MIS, ADp, CDp, Time, G, GV, US, param_file, diag
   real :: wave_speed_tol      ! The fractional tolerance for finding the wave speeds [nondim]
   logical :: better_speed_est ! If true, use a more robust estimate of the first
                               ! mode wave speed as the starting point for iterations.
+  logical :: split            ! True if using the barotropic-baroclinic split algorithm
   logical :: use_temperature, adiabatic
   logical :: default_2018_answers, remap_answers_2018
   integer :: isd, ied, jsd, jed, IsdB, IedB, JsdB, JedB, nz, nkml, nkbl
@@ -1568,6 +1589,7 @@ subroutine MOM_diagnostics_init(MIS, ADp, CDp, Time, G, GV, US, param_file, diag
                  "If true, use the order of arithmetic and expressions that recover the "//&
                  "answers from the end of 2018.  Otherwise, use updated and more robust "//&
                  "forms of the same expressions.", default=default_2018_answers)
+  call get_param(param_file, mdl, "SPLIT", split, default=.true., do_not_log=.true.)
 
   if (GV%Boussinesq) then
     thickness_units = "m" ; flux_units = "m3 s-1" ; convert_H = GV%H_to_m
@@ -1779,6 +1801,13 @@ subroutine MOM_diagnostics_init(MIS, ADp, CDp, Time, G, GV, US, param_file, diag
       'm3 s-3', conversion=GV%H_to_m*(US%L_T_to_m_s**2)*US%s_to_T)
   if (CS%id_PE_to_KE>0) call safe_alloc_ptr(CS%PE_to_KE,isd,ied,jsd,jed,nz)
 
+  if (split) then
+    CS%id_KE_BT = register_diag_field('ocean_model', 'KE_BT', diag%axesTL, Time, &
+        'Barotropic contribution to Kinetic Energy', &
+        'm3 s-3', conversion=GV%H_to_m*(US%L_T_to_m_s**2)*US%s_to_T)
+    if (CS%id_KE_BT>0) call safe_alloc_ptr(CS%KE_BT,isd,ied,jsd,jed,nz)
+  endif
+
   CS%id_KE_Coradv = register_diag_field('ocean_model', 'KE_Coradv', diag%axesTL, Time, &
       'Kinetic Energy Source from Coriolis and Advection', &
       'm3 s-3', conversion=GV%H_to_m*(US%L_T_to_m_s**2)*US%s_to_T)
@@ -2192,9 +2221,9 @@ subroutine set_dependent_diagnostics(MIS, ADp, CDp, G, CS)
   IsdB = G%IsdB ; IedB = G%IedB ; JsdB = G%JsdB ; JedB = G%JedB
 
   if (associated(CS%dKE_dt) .or. associated(CS%PE_to_KE) .or. &
-      associated(CS%KE_CorAdv) .or. associated(CS%KE_adv) .or. &
-      associated(CS%KE_visc) .or. associated(CS%KE_horvisc) .or. &
-      associated(CS%KE_dia)) &
+      associated(CS%KE_BT) .or. associated(CS%KE_CorAdv) .or. &
+      associated(CS%KE_adv) .or. associated(CS%KE_visc) .or. &
+      associated(CS%KE_horvisc) .or. associated(CS%KE_dia)) &
     call safe_alloc_ptr(CS%KE,isd,ied,jsd,jed,nz)
 
   if (associated(CS%dKE_dt)) then
@@ -2246,6 +2275,7 @@ subroutine MOM_diagnostics_end(CS, ADp)
   if (associated(CS%KE))         deallocate(CS%KE)
   if (associated(CS%dKE_dt))     deallocate(CS%dKE_dt)
   if (associated(CS%PE_to_KE))   deallocate(CS%PE_to_KE)
+  if (associated(CS%KE_BT))      deallocate(CS%KE_BT)
   if (associated(CS%KE_Coradv))  deallocate(CS%KE_Coradv)
   if (associated(CS%KE_adv))     deallocate(CS%KE_adv)
   if (associated(CS%KE_visc))    deallocate(CS%KE_visc)

src/parameterizations/lateral/MOM_hor_visc.F90

@@ -193,6 +193,7 @@ module MOM_hor_visc
   integer :: id_sh_xy_q = -1,    id_sh_xx_h      = -1
   integer :: id_FrictWork = -1, id_FrictWorkIntz = -1
   integer :: id_FrictWork_GME = -1
+  integer :: id_normstress = -1, id_shearstress = -1
   !>@}
 
 
@@ -307,8 +308,8 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US,
     vort_xy_q, & ! vertical vorticity at corner points [T-1 ~> s-1]
     sh_xy_q,   & ! horizontal shearing strain at corner points [T-1 ~> s-1]
     GME_coeff_q, &  !< GME coeff. at q-points [L2 T-1 ~> m2 s-1]
-    max_diss_rate_q ! maximum possible energy dissipated by lateral friction [L2 T-3 ~> m2 s-3]
-
+    max_diss_rate_q, & ! maximum possible energy dissipated by lateral friction [L2 T-3 ~> m2 s-3]
+    ShSt         ! A diagnostic array of shear stress [T-1 ~> s-1].
   real, dimension(SZIB_(G),SZJ_(G),SZK_(G)+1) :: &
     KH_u_GME  !< interface height diffusivities in u-columns [L2 T-1 ~> m2 s-1]
   real, dimension(SZI_(G),SZJB_(G),SZK_(G)+1) :: &
@@ -320,7 +321,8 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US,
     FrictWork, &     ! work done by MKE dissipation mechanisms [R L2 T-3 ~> W m-2]
     FrictWork_GME, & ! work done by GME [R L2 T-3 ~> W m-2]
     div_xx_h,      & ! horizontal divergence [T-1 ~> s-1]
-    sh_xx_h          ! horizontal tension (du/dx - dv/dy) including metric terms [T-1 ~> s-1]
+    sh_xx_h,       & ! horizontal tension (du/dx - dv/dy) including metric terms [T-1 ~> s-1]
+    NoSt             ! A diagnostic array of normal stress [T-1 ~> s-1].
   real, dimension(SZI_(G),SZJ_(G),SZK_(G)) :: &
     grid_Re_Kh, &    !< Grid Reynolds number for Laplacian horizontal viscosity at h points [nondim]
     grid_Re_Ah, &    !< Grid Reynolds number for Biharmonic horizontal viscosity at h points [nondim]
@@ -494,7 +496,7 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US,
   !$OMP   diffu, diffv, max_diss_rate_h, max_diss_rate_q, &
   !$OMP   Kh_h, Kh_q, Ah_h, Ah_q, FrictWork, FrictWork_GME, &
   !$OMP   div_xx_h, sh_xx_h, vort_xy_q, sh_xy_q, GME_coeff_h, GME_coeff_q, &
-  !$OMP   TD, KH_u_GME, KH_v_GME, grid_Re_Kh, grid_Re_Ah &
+  !$OMP   TD, KH_u_GME, KH_v_GME, grid_Re_Kh, grid_Re_Ah, NoSt, ShSt &
   !$OMP ) &
   !$OMP private( &
   !$OMP   i, j, k, n, &
@@ -524,6 +526,7 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US,
       dvdy(i,j) = CS%DX_dyT(i,j)*(G%IdxCv(i,J) * v(i,J,k) - &
                                   G%IdxCv(i,J-1) * v(i,J-1,k))
       sh_xx(i,j) = dudx(i,j) - dvdy(i,j)
+      if (CS%id_normstress > 0) NoSt(i,j,k) = sh_xx(i,j)
     enddo ; enddo
 
     ! Components for the shearing strain
@@ -671,10 +674,12 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US,
     if (CS%no_slip) then
       do J=js-2,Jeq+1 ; do I=is-2,Ieq+1
         sh_xy(I,J) = (2.0-G%mask2dBu(I,J)) * ( dvdx(I,J) + dudy(I,J) )
+        if (CS%id_shearstress > 0) ShSt(I,J,k) = sh_xy(I,J)
       enddo ; enddo
     else
       do J=js-2,Jeq+1 ; do I=is-2,Ieq+1
         sh_xy(I,J) = G%mask2dBu(I,J) * ( dvdx(I,J) + dudy(I,J) )
+        if (CS%id_shearstress > 0) ShSt(I,J,k) = sh_xy(I,J)
       enddo ; enddo
     endif
 
@@ -1338,6 +1343,8 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US,
   enddo ! end of k loop
 
   ! Offer fields for diagnostic averaging.
+  if (CS%id_normstress > 0) call post_data(CS%id_normstress, NoSt, CS%diag)
+  if (CS%id_shearstress > 0) call post_data(CS%id_shearstress, ShSt, CS%diag)
   if (CS%id_diffu>0)     call post_data(CS%id_diffu, diffu, CS%diag)
   if (CS%id_diffv>0)     call post_data(CS%id_diffv, diffv, CS%diag)
   if (CS%id_FrictWork>0) call post_data(CS%id_FrictWork, FrictWork, CS%diag)
@@ -2074,6 +2081,10 @@ subroutine hor_visc_init(Time, G, US, param_file, diag, CS, MEKE, ADp)
     endif
   endif
   ! Register fields for output from this module.
+  CS%id_normstress = register_diag_field('ocean_model', 'NoSt', diag%axesTL, Time, &
+      'Normal Stress', 's-1', conversion=US%s_to_T)
+  CS%id_shearstress = register_diag_field('ocean_model', 'ShSt', diag%axesBL, Time, &
+      'Shear Stress', 's-1', conversion=US%s_to_T)
   CS%id_diffu = register_diag_field('ocean_model', 'diffu', diag%axesCuL, Time, &
       'Zonal Acceleration from Horizontal Viscosity', 'm s-2', conversion=US%L_T2_to_m_s2)
   CS%id_diffv = register_diag_field('ocean_model', 'diffv', diag%axesCvL, Time, &

src/user/basin_builder.F90

@@ -94,6 +94,18 @@ subroutine basin_builder_topography(D, G, param_file, max_depth)
         lat = G%geoLatT(i,j)
         D(i,j) = min( D(i,j), NS_scurve_ridge(lon, lat, pars(1), pars(2), pars(3), pars(4), pars(5)) )
       enddo ; enddo
+    elseif (trim(lowercase(funcs)) == 'angled_coast') then
+      call get_param(param_file, mdl, pname2, pars(1:4), &
+                     "ANGLED_COAST parameters: longitude intersection with Equator, "//&
+                     "latitude intersection with Prime Meridian, footprint radius, shelf depth.", &
+                     units="degrees_E,degrees_N,degrees,m", &
+                     fail_if_missing=.true.)
+      pars(4) = pars(4) / max_depth
+      do j=G%jsc,G%jec ; do i=G%isc,G%iec
+        lon = G%geoLonT(i,j)
+        lat = G%geoLatT(i,j)
+        D(i,j) = min( D(i,j), angled_coast(lon, lat, pars(1), pars(2), pars(3), pars(4)) )
+      enddo ; enddo
     elseif (trim(lowercase(funcs)) == 'ew_coast') then
       call get_param(param_file, mdl, pname2, pars(1:5), &
                      "EW_COAST parameters: latitude, starting longitude, "//&
@@ -211,6 +223,21 @@ real function dist_line_fixed_y(x, y, x0, x1, y0)
   dist_line_fixed_y = dist_line_fixed_x(y, x, y0, x0, x1)
 end function dist_line_fixed_y
 
+!> An "angled coast profile".
+real function angled_coast(lon, lat, lon_eq, lat_mer, dr, sh)
+  real, intent(in) :: lon     !< Longitude [degrees_E]
+  real, intent(in) :: lat     !< Latitude [degrees_N]
+  real, intent(in) :: lon_eq  !< Longitude intersection with Equator [degrees_E]
+  real, intent(in) :: lat_mer !< Latitude intersection with Prime Meridian [degrees_N]
+  real, intent(in) :: dr      !< "Radius" of coast profile [degrees]
+  real, intent(in) :: sh      !< depth of shelf as fraction of full depth [nondim]
+  real :: r
+
+  r = 1/sqrt( lat_mer*lat_mer + lon_eq*lon_eq )
+  r = r * ( lat_mer*lon + lon_eq*lat - lon_eq*lat_mer)
+  angled_coast = cstprof(r, 0., dr, 0.125, 0.125, 0.5, sh)
+end function angled_coast
+
 !> A "coast profile" applied in an N-S line from lonC,lat0 to lonC,lat1.
 real function NS_coast(lon, lat, lonC, lat0, lat1, dlon, sh)
   real, intent(in) :: lon     !< Longitude [degrees_E]