Merge branch 'dev/gfdl' into tc-reconfig

.testing/Makefile

@@ -1,8 +1,12 @@
 SHELL = bash
-MPIRUN ?= mpirun
 
+# User-defined configuration
 -include config.mk
 
+# Default configurations
+MPIRUN ?= mpirun
+DO_REPRO_TESTS ?= true
+
 #---
 # Dependencies
 BASE = $(dir $(abspath $(lastword $(MAKEFILE_LIST))))/..
@@ -36,38 +40,47 @@ MKMF_TEMPLATE ?= $(DEPS)/mkmf/templates/ncrc-gnu.mk
 # Test configuration
 
 # Executables
-BUILDS = symmetric asymmetric repro
+BUILDS = symmetric asymmetric repro openmp
 CONFIGS := $(wildcard tc*)
-TESTS = grids layouts restarts repros nans dims
+TESTS = grids layouts restarts nans dims openmps
+
+# REPRO tests enable reproducibility with optimization, and often do not match
+# the DEBUG results in older GCCs and vendor compilers, so we can optionally
+# disable them.
+ifeq ($(DO_REPRO_TESTS), true)
+  BUILDS += repro
+  TESTS += repros
+endif
 
 # The following variables are configured by Travis:
 #   DO_REGRESSION_TESTS: true if $(TRAVIS_PULL_REQUEST) is a PR number
 #   MOM_TARGET_SLUG: TRAVIS_REPO_SLUG
 #   MOM_TARGET_LOCAL_BRANCH: TRAVIS_BRANCH
-#
-# These are set to true by Travis if testing a pull request
+
+# These are set to true by our Travis configuration if testing a pull request
 DO_REGRESSION_TESTS ?=
 REPORT_COVERAGE ?=
 
 ifeq ($(DO_REGRESSION_TESTS), true)
-	BUILDS += target
-	TEST += regressions
+  BUILDS += target
+  TESTS += regressions
 
-	MOM_TARGET_SLUG ?= NOAA-GFDL/MOM6
-	MOM_TARGET_URL ?= https://github.com/$(MOM_TARGET_SLUG)
+  MOM_TARGET_SLUG ?= NOAA-GFDL/MOM6
+  MOM_TARGET_URL ?= https://github.com/$(MOM_TARGET_SLUG)
 
-	MOM_TARGET_LOCAL_BRANCH ?= dev/gfdl
-	MOM_TARGET_BRANCH := origin/$(MOM_TARGET_LOCAL_BRANCH)
+  MOM_TARGET_LOCAL_BRANCH ?= dev/gfdl
+  MOM_TARGET_BRANCH := origin/$(MOM_TARGET_LOCAL_BRANCH)
 
-	TARGET_CODEBASE = $(BUILD)/target_codebase
+  TARGET_CODEBASE = $(BUILD)/target_codebase
 else
-	MOM_TARGET_URL =
-	MOM_TARGET_BRANCH =
-	TARGET_CODEBASE =
+  MOM_TARGET_URL =
+  MOM_TARGET_BRANCH =
+  TARGET_CODEBASE =
 endif
 
 SOURCE = $(wildcard $(BASE)/src/*/*.F90 $(BASE)/src/*/*/*.F90 $(BASE)/config_src/solo_driver/*.F90)
 
+
 #---
 # Rules
 
@@ -85,6 +98,7 @@ $(BUILD)/target/MOM6: MOMFLAGS=NETCDF=3 DEBUG=1
 $(BUILD)/symmetric/MOM6: MOMFLAGS=NETCDF=3 DEBUG=1 $(COVFLAG)
 $(BUILD)/asymmetric/MOM6: MOMFLAGS=NETCDF=3 DEBUG=1
 $(BUILD)/repro/MOM6: MOMFLAGS=NETCDF=3 REPRO=1
+$(BUILD)/openmp/MOM6: MOMFLAGS=NETCDF=3 OPENMP=1
 
 $(BUILD)/asymmetric/path_names: GRID_SRC=config_src/dynamic
 $(BUILD)/%/path_names: GRID_SRC=config_src/dynamic_symmetric
@@ -170,44 +184,48 @@ test.grids: $(foreach c,$(filter-out tc3,$(CONFIGS)),$(c).grid $(c).grid.diag)
 test.layouts: $(foreach c,$(CONFIGS),$(c).layout $(c).layout.diag)
 test.restarts: $(foreach c,$(CONFIGS),$(c).restart)
 test.repros: $(foreach c,$(CONFIGS),$(c).repro $(c).repro.diag)
+test.openmps: $(foreach c,$(CONFIGS),$(c).openmp $(c).openmp.diag)
 test.nans: $(foreach c,$(CONFIGS),$(c).nan $(c).nan.diag)
 test.dims: $(foreach c,$(CONFIGS),$(foreach d,t l h z,$(c).dim.$(d) $(c).dim.$(d).diag))
 
-# NOTE: chksum_diag return code of cmp is currently ignored since many fail!
 define CMP_RULE
 .PRECIOUS: $(foreach b,$(2),results/%/ocean.stats.$(b))
 %.$(1): $(foreach b,$(2),results/%/ocean.stats.$(b))
-	cmp $$^
+	cmp $$^ || diff $$^
 
 .PRECIOUS: $(foreach b,$(2),results/%/chksum_diag.$(b))
 %.$(1).diag: $(foreach b,$(2),results/%/chksum_diag.$(b))
-	cmp $$^ || true
+	cmp $$^ || diff $$^
 endef
 
 $(eval $(call CMP_RULE,regression,symmetric target))
 $(eval $(call CMP_RULE,grid,symmetric asymmetric))
 $(eval $(call CMP_RULE,layout,symmetric layout))
 $(eval $(call CMP_RULE,repro,symmetric repro))
+$(eval $(call CMP_RULE,openmp,symmetric openmp))
 $(eval $(call CMP_RULE,nan,symmetric nan))
 $(foreach d,t l h z,$(eval $(call CMP_RULE,dim.$(d),symmetric dim.$(d))))
 
 # Restart tests only compare the final stat record
 .PRECIOUS: $(foreach b,symmetric restart,results/%/ocean.stats.$(b))
 %.restart: $(foreach b,symmetric restart,results/%/ocean.stats.$(b))
-	cmp $(foreach f,$^,<(tr -s ' ' < $(f) | cut -d ' ' -f3- | tail -n 1))
+	cmp $(foreach f,$^,<(tr -s ' ' < $(f) | cut -d ' ' -f3- | tail -n 1)) \
+		|| diff $^
 
 # TODO: chksum_diag parsing of restart files
 
 
 #---
 # Test run output files
 
-# Simple function for generalized Slurm (srun) and OpenMPI (mpirun) support
+# Generalized MPI environment variable support
 # $(1): Environment variables
-ifeq ($(MPIRUN), srun)
-MPIRUN_CMD=$(1) $(MPIRUN)
+ifeq ($(shell $(MPIRUN) -x tmp=1 true 2> /dev/null ; echo $$?), 0)
+  MPIRUN_CMD=$(MPIRUN) $(if $(1),-x $(1),)
+else ifeq ($(shell $(MPIRUN) -env tmp=1 true 2> /dev/null ; echo $$?), 0)
+  MPIRUN_CMD=$(MPIRUN) $(if $(1),-env $(1),)
 else
-MPIRUN_CMD=$(MPIRUN) $(if $(1),-x $(1),)
+  MPIRUN_CMD=$(1) $(MPIRUN)
 endif
 
 # Rule to build results/<tc>/{ocean.stats,chksum_diag}.<tag>
@@ -243,6 +261,7 @@ $(eval $(call STAT_RULE,symmetric,symmetric,$(REPORT_COVERAGE),,,1))
 $(eval $(call STAT_RULE,asymmetric,asymmetric,,,,1))
 $(eval $(call STAT_RULE,target,target,,,,1))
 $(eval $(call STAT_RULE,repro,repro,,,,1))
+$(eval $(call STAT_RULE,openmp,openmp,,,,1))
 $(eval $(call STAT_RULE,layout,symmetric,,LAYOUT=2$(,)1,,2))
 $(eval $(call STAT_RULE,nan,symmetric,,,MALLOC_PERTURB_=256,1))
 $(eval $(call STAT_RULE,dim.t,symmetric,,T_RESCALE_POWER=11,,1))
@@ -255,7 +274,7 @@ results/%/ocean.stats.restart: ../build/symmetric/MOM6
 	rm -rf work/$*/restart
 	mkdir -p work/$*/restart
 	cp -rL $*/* work/$*/restart
-	mkdir work/$*/restart/RESTART
+	mkdir -p work/$*/restart/RESTART
 	# Generate the half-period input namelist
 	# TODO: Assumes runtime set by DAYMAX, will fail if set by input.nml
 	cd work/$*/restart \

src/ALE/MOM_ALE.F90

@@ -799,12 +799,12 @@ subroutine remap_all_state_vars(CS_remapping, CS_ALE, G, GV, h_old, h_new, Reg,
 
         ! Intermediate steps for tendency of tracer concentration and tracer content.
         if (present(dt)) then
-          if (Tr%id_remap_conc>0) then
+          if (Tr%id_remap_conc > 0) then
             do k=1,GV%ke
-              work_conc(i,j,k) = (u_column(k)    - Tr%t(i,j,k)      ) * Idt
+              work_conc(i,j,k) = (u_column(k) - Tr%t(i,j,k)) * Idt
             enddo
           endif
-          if (Tr%id_remap_cont>0. .or. Tr%id_remap_cont_2d>0) then
+          if (Tr%id_remap_cont > 0 .or. Tr%id_remap_cont_2d > 0) then
             do k=1,GV%ke
               work_cont(i,j,k) = (u_column(k)*h2(k) - Tr%t(i,j,k)*h1(k)) * Idt
             enddo

src/core/MOM.F90

@@ -1674,8 +1674,8 @@ subroutine initialize_MOM(Time, Time_init, param_file, dirs, CS, restart_CSp, &
                  "faster by eliminating subroutine calls.", default=.false.)
   call get_param(param_file, "MOM", "DO_DYNAMICS", CS%do_dynamics, &
                  "If False, skips the dynamics calls that update u & v, as well as "//&
-                 "the gravity wave adjustment to h. This is a fragile feature and "//&
-                 "thus undocumented.", default=.true., do_not_log=.true. )
+                 "the gravity wave adjustment to h. This may be a fragile feature, "//&
+                 "but can be useful during development", default=.true.)
   call get_param(param_file, "MOM", "ADVECT_TS", advect_TS, &
                  "If True, advect temperature and salinity horizontally "//&
                  "If False, T/S are registered for advection. "//&

src/core/MOM_PressureForce_blocked_AFV.F90

@@ -327,7 +327,7 @@ subroutine PressureForce_blk_AFV_nonBouss(h, tv, PFu, PFv, G, GV, US, CS, p_atm,
   ! geopotentials will not now be linear at the sub-grid-scale.  Doing this
   ! ensures no motion with flat isopycnals, even with a nonlinear equation of state.
 !$OMP parallel do default(none) shared(nz,za,G,GV,dza,intx_dza,h,PFu, &
-!$OMP                                  intp_dza,p,dp_neglect,inty_dza,PFv,CS,dM) &
+!$OMP                                  intp_dza,p,dp_neglect,inty_dza,PFv,CS,dM,US) &
 !$OMP                          private(is_bk,ie_bk,js_bk,je_bk,Isq_bk,Ieq_bk,Jsq_bk, &
 !$OMP                                  Jeq_bk,ioff_bk,joff_bk,i,j,za_bk,intx_za_bk,  &
 !$OMP                                  inty_za_bk,dp_bk)

src/core/MOM_barotropic.F90

@@ -1262,9 +1262,9 @@ subroutine btstep(U_in, V_in, eta_in, dt, bc_accel_u, bc_accel_v, forces, pbce,
 !$OMP                               find_etaav,jsvf,jevf,isvf,ievf,eta_sum,eta_wtd, &
 !$OMP                               ubt_sum,uhbt_sum,PFu_bt_sum,Coru_bt_sum,ubt_wtd,&
 !$OMP                               ubt_trans,vbt_sum,vhbt_sum,PFv_bt_sum,          &
-!$OMP                               Corv_bt_sum,vbt_wtd,vbt_trans,eta_src,dt,dtbt,  &
+!$OMP                               Corv_bt_sum,vbt_wtd,vbt_trans,eta_src,dt_in_T,dtbt,  &
 !$OMP                               Rayleigh_u, Rayleigh_v,                         &
-!$OMP                               use_BT_Cont,BTCL_u,uhbt0,BTCL_v,vhbt0,eta,Idt)  &
+!$OMP                               use_BT_Cont,BTCL_u,uhbt0,BTCL_v,vhbt0,eta,Idt,US)  &
 !$OMP                       private(u_max_cor,v_max_cor,eta_cor_max,Htot)
   !$OMP do
   do j=js-1,je+1 ; do I=is-1,ie ; av_rem_u(I,j) = 0.0 ; enddo ; enddo

src/core/MOM_continuity_PPM.F90

@@ -300,7 +300,7 @@ subroutine zonal_mass_flux(u, h_in, uh, dt_in_T, G, GV, US, CS, LB, uhbt, OBC, &
 
   call cpu_clock_begin(id_clock_correct)
 !$OMP parallel do default(none) shared(ish,ieh,jsh,jeh,nz,u,h_in,h_L,h_R,use_visc_rem,visc_rem_u,  &
-!$OMP                                  uh,dt,G,GV,CS,local_specified_BC,OBC,uhbt,set_BT_cont,    &
+!$OMP                                  uh,dt_in_T,US,G,GV,CS,local_specified_BC,OBC,uhbt,set_BT_cont,    &
 !$OMP                                  CFL_dt,I_dt,u_cor,BT_cont, local_Flather_OBC) &
 !$OMP                          private(do_I,duhdu,du,du_max_CFL,du_min_CFL,uh_tot_0,duhdu_tot_0, &
 !$OMP                                  is_simple,FAuI,visc_rem_max,I_vrm,du_lim,dx_E,dx_W,any_simple_OBC ) &
@@ -1099,7 +1099,7 @@ subroutine meridional_mass_flux(v, h_in, vh, dt_in_T, G, GV, US, CS, LB, vhbt, O
 
   call cpu_clock_begin(id_clock_correct)
 !$OMP parallel do default(none) shared(ish,ieh,jsh,jeh,nz,v,h_in,h_L,h_R,vh,use_visc_rem, &
-!$OMP                                  visc_rem_v,dt,G,GV,CS,local_specified_BC,OBC,vhbt, &
+!$OMP                                  visc_rem_v,dt_in_T,US,G,GV,CS,local_specified_BC,OBC,vhbt, &
 !$OMP                                  set_BT_cont,CFL_dt,I_dt,v_cor,BT_cont, local_Flather_OBC ) &
 !$OMP                          private(do_I,dvhdv,dv,dv_max_CFL,dv_min_CFL,vh_tot_0,    &
 !$OMP                                  dvhdv_tot_0,visc_rem_max,I_vrm,dv_lim,dy_N,      &

src/diagnostics/MOM_diagnostics.F90

@@ -623,7 +623,7 @@ subroutine calculate_diagnostic_fields(u, v, h, uh, vh, tv, ADp, CDp, p_surf, &
     call wave_speed(h, tv, G, GV, US, CS%cg1, CS%wave_speed_CSp)
     if (CS%id_cg1>0) call post_data(CS%id_cg1, CS%cg1, CS%diag)
     if (CS%id_Rd1>0) then
-!$OMP parallel do default(none) shared(is,ie,js,je,G,CS) &
+!$OMP parallel do default(none) shared(is,ie,js,je,G,CS,US) &
 !$OMP                          private(f2_h,mag_beta)
       do j=js,je ; do i=is,ie
         ! Blend the equatorial deformation radius with the standard one.
@@ -672,7 +672,7 @@ subroutine calculate_diagnostic_fields(u, v, h, uh, vh, tv, ADp, CDp, p_surf, &
     endif
     if (CS%id_cg_ebt>0) call post_data(CS%id_cg_ebt, CS%cg1, CS%diag)
     if (CS%id_Rd_ebt>0) then
-!$OMP parallel do default(none) shared(is,ie,js,je,G,CS) &
+!$OMP parallel do default(none) shared(is,ie,js,je,G,CS,US) &
 !$OMP                          private(f2_h,mag_beta)
       do j=js,je ; do i=is,ie
         ! Blend the equatorial deformation radius with the standard one.

src/equation_of_state/MOM_EOS.F90

@@ -82,6 +82,11 @@ module MOM_EOS
   module procedure calculate_TFreeze_scalar, calculate_TFreeze_array
 end interface calculate_TFreeze
 
+!> Calculates the compressibility of water from T, S, and P
+interface calculate_compress
+  module procedure calculate_compress_scalar, calculate_compress_array
+end interface calculate_compress
+
 !> A control structure for the equation of state
 type, public :: EOS_type ; private
   integer :: form_of_EOS = 0 !< The equation of state to use.
@@ -523,16 +528,16 @@ subroutine calculate_specific_vol_derivs(T, S, pressure, dSV_dT, dSV_dS, start,
 end subroutine calculate_specific_vol_derivs
 
 !> Calls the appropriate subroutine to calculate the density and compressibility for 1-D array inputs.
-subroutine calculate_compress(T, S, pressure, rho, drho_dp, start, npts, EOS)
-  real, dimension(:), intent(in)  :: T !< Potential temperature referenced to the surface [degC]
-  real, dimension(:), intent(in)  :: S !< Salinity [ppt]
+subroutine calculate_compress_array(T, S, pressure, rho, drho_dp, start, npts, EOS)
+  real, dimension(:), intent(in)  :: T        !< Potential temperature referenced to the surface [degC]
+  real, dimension(:), intent(in)  :: S        !< Salinity [PSU]
   real, dimension(:), intent(in)  :: pressure !< Pressure [Pa]
-  real, dimension(:), intent(out) :: rho !< In situ density [kg m-3].
-  real, dimension(:), intent(out) :: drho_dp !< The partial derivative of density with pressure
-                                     !! (also the inverse of the square of sound speed) in s2 m-2.
-  integer,            intent(in)  :: start !< Starting index within the array
-  integer,            intent(in)  :: npts !< The number of values to calculate
-  type(EOS_type),     pointer     :: EOS !< Equation of state structure
+  real, dimension(:), intent(out) :: rho      !< In situ density [kg m-3].
+  real, dimension(:), intent(out) :: drho_dp   !< The partial derivative of density with pressure
+                                     !! (also the inverse of the square of sound speed) [s2 m-2].
+  integer,            intent(in)  :: start    !< Starting index within the array
+  integer,            intent(in)  :: npts     !< The number of values to calculate
+  type(EOS_type),     pointer     :: EOS      !< Equation of state structure
 
   if (.not.associated(EOS)) call MOM_error(FATAL, &
     "calculate_compress called with an unassociated EOS_type EOS.")
@@ -554,8 +559,29 @@ subroutine calculate_compress(T, S, pressure, rho, drho_dp, start, npts, EOS)
            "calculate_compress: EOS%form_of_EOS is not valid.")
   end select
 
-end subroutine calculate_compress
+end subroutine calculate_compress_array
+
+!> Calculate density and compressibility for a scalar. This just promotes the scalar to an array with a singleton
+!! dimension and calls calculate_compress_array
+subroutine calculate_compress_scalar(T, S, pressure, rho, drho_dp, EOS)
+  real, intent(in)        :: T        !< Potential temperature referenced to the surface (degC)
+  real, intent(in)        :: S        !< Salinity (PSU)
+  real, intent(in)        :: pressure !< Pressure (Pa)
+  real, intent(out)       :: rho      !< In situ density in kg m-3.
+  real, intent(out)       :: drho_dp  !< The partial derivative of density with pressure
+                                      !! (also the inverse of the square of sound speed) in s2 m-2.
+  type(EOS_type), pointer :: EOS      !< Equation of state structure
+
+  real, dimension(1) :: Ta, Sa, pa, rhoa, drho_dpa
+
+  if (.not.associated(EOS)) call MOM_error(FATAL, &
+    "calculate_compress called with an unassociated EOS_type EOS.")
+  Ta(1) = T ; Sa(1) = S; pa(1) = pressure
+
+  call calculate_compress_array(Ta, Sa, pa, rhoa, drho_dpa, 1, 1, EOS)
+  rho = rhoa(1) ; drho_dp = drho_dpa(1)
 
+end subroutine calculate_compress_scalar
 !> Calls the appropriate subroutine to alculate analytical and nearly-analytical
 !! integrals in pressure across layers of geopotential anomalies, which are
 !! required for calculating the finite-volume form pressure accelerations in a

src/framework/MOM_checksums.F90

@@ -1822,8 +1822,11 @@ subroutine chk_sum_msg3(fmsg, aMean, aMin, aMax, mesg, iounit)
   real,             intent(in) :: aMax !< The maximum value of the array
   integer,          intent(in) :: iounit !< Checksum logger IO unit
 
+  ! NOTE: We add zero to aMin and aMax to remove any negative zeros.
+  ! This is due to inconsistencies of signed zero in local vs MPI calculations.
+
   if (is_root_pe()) write(iounit, '(A,3(A,ES25.16,1X),A)') &
-    fmsg, " mean=", aMean, "min=", aMin, "max=", aMax, trim(mesg)
+    fmsg, " mean=", aMean, "min=", (0. + aMin), "max=", (0. + aMax), trim(mesg)
 end subroutine chk_sum_msg3
 
 !> MOM_checksums_init initializes the MOM_checksums module. As it happens, the
@@ -1849,20 +1852,12 @@ end subroutine chksum_error
 !> Does a bitcount of a number by first casting to an integer and then using BTEST
 !! to check bit by bit
 integer function bitcount(x)
-  real :: x !< Number to be bitcount
-
-  ! Local variables
-  integer(kind(x)) :: y !< Store the integer representation of the memory used by x
-  integer :: bit
+  real, intent(in) :: x !< Number to be bitcount
 
-  bitcount = 0
-  y = transfer(x,y)
-
-  ! Fortran standard says that bit indexing start at 0
-  do bit = 0, bit_size(y)-1
-    if (BTEST(y,bit)) bitcount = bitcount+1
-  enddo
+  integer, parameter :: xk = kind(x)  !< Kind type of x
 
+  ! NOTE: Assumes that reals and integers of kind=xk are the same size
+  bitcount = popcnt(transfer(x, 1_xk))
 end function bitcount
 
 end module MOM_checksums

src/parameterizations/lateral/MOM_MEKE.F90

@@ -128,7 +128,8 @@ subroutine step_forward_MEKE(MEKE, h, SN_u, SN_v, visc, dt, G, GV, US, CS, hu, h
     del4MEKE, &     ! Time-integrated MEKE tendency arising from the biharmonic of MEKE [L2 T-2 ~> m2 s-2].
     LmixScale, &    ! Eddy mixing length [L ~> m].
     barotrFac2, &   ! Ratio of EKE_barotropic / EKE [nondim]
-    bottomFac2      ! Ratio of EKE_bottom / EKE [nondim]
+    bottomFac2, &   ! Ratio of EKE_bottom / EKE [nondim]
+    tmp             ! Temporary variable for diagnostic computation
 
   real, dimension(SZIB_(G),SZJ_(G)) :: &
     MEKE_uflux, &   ! The zonal advective and diffusive flux of MEKE with different units in different
@@ -575,10 +576,29 @@ subroutine step_forward_MEKE(MEKE, h, SN_u, SN_v, visc, dt, G, GV, US, CS, hu, h
     endif
 
     ! Offer fields for averaging.
+
+    if (any([CS%id_Ue, CS%id_Ub, CS%id_Ut] > 0)) &
+      tmp(:,:) = 0.
+
     if (CS%id_MEKE>0) call post_data(CS%id_MEKE, MEKE%MEKE, CS%diag)
-    if (CS%id_Ue>0) call post_data(CS%id_Ue, sqrt(max(0.,2.0*MEKE%MEKE)), CS%diag)
-    if (CS%id_Ub>0) call post_data(CS%id_Ub, sqrt(max(0.,2.0*MEKE%MEKE*bottomFac2)), CS%diag)
-    if (CS%id_Ut>0) call post_data(CS%id_Ut, sqrt(max(0.,2.0*MEKE%MEKE*barotrFac2)), CS%diag)
+    if (CS%id_Ue>0) then
+      do j=js,je ; do i=is,ie
+        tmp(i,j) = sqrt(max(0., 2. * MEKE%MEKE(i,j)))
+      enddo ; enddo
+      call post_data(CS%id_Ue, tmp, CS%diag)
+    endif
+    if (CS%id_Ub>0) then
+      do j=js,je ; do i=is,ie
+        tmp(i,j) = sqrt(max(0., 2. * MEKE%MEKE(i,j) * bottomFac2(i,j)))
+      enddo ; enddo
+      call post_data(CS%id_Ub, tmp, CS%diag)
+    endif
+    if (CS%id_Ut>0) then
+      do j=js,je ; do i=is,ie
+        tmp(i,j) = sqrt(max(0., 2. * MEKE%MEKE(i,j) * barotrFac2(i,j)))
+      enddo ; enddo
+      call post_data(CS%id_Ut, tmp, CS%diag)
+    endif
     if (CS%id_Kh>0) call post_data(CS%id_Kh, MEKE%Kh, CS%diag)
     if (CS%id_Ku>0) call post_data(CS%id_Ku, MEKE%Ku, CS%diag)
     if (CS%id_Au>0) call post_data(CS%id_Au, MEKE%Au, CS%diag)