Merge branch 'dev/ncar' into dev-master-candidate-2018-06-03

config_src/mct_driver/ocn_comp_mct.F90

@@ -1774,7 +1774,7 @@ subroutine update_ocean_model(OS, Ocean_sfc, time_start_update, &
   call set_net_mass_forcing(OS%fluxes, OS%forces, OS%grid)
 
   if (OS%nstep==0) then
-    call finish_MOM_initialization(OS%Time, OS%dirs, OS%MOM_CSp, S%restart_CSp)
+    call finish_MOM_initialization(OS%Time, OS%dirs, OS%MOM_CSp, OS%restart_CSp)
   endif
 
   call disable_averaging(OS%diag)

src/parameterizations/vertical/MOM_KPP.F90

@@ -88,7 +88,8 @@ module MOM_KPP
   character(len=30) :: MatchTechnique  !< Method used in CVMix for setting diffusivity and NLT profile functions
   integer :: NLT_shape                 !< MOM6 over-ride of CVMix NLT shape function
   logical :: applyNonLocalTrans        !< If True, apply non-local transport to heat and scalars
-  logical :: smoothBLD                 !< If True, apply a 1-1-4-1-1 Laplacian filter one time on HBLT.
+  integer :: n_smooth                  !< Number of times smoothing operator is applied on OBLdepth.
+  logical :: deepen_only               !< If true, apply OBLdepth smoothing at a cell only if the OBLdepth gets deeper.
   logical :: KPPzeroDiffusivity        !< If True, will set diffusivity and viscosity from KPP to zero
                                        !! for testing purposes.
   logical :: KPPisAdditive             !< If True, will add KPP diffusivity to initial diffusivity.
@@ -134,10 +135,12 @@ module MOM_KPP
   integer :: id_NLT_temp_budget = -1
   integer :: id_NLT_saln_budget = -1
   integer :: id_EnhK = -1, id_EnhW = -1, id_EnhVt2 = -1
+  integer :: id_OBLdepth_original = -1
 
 
   ! Diagnostics arrays
   real, allocatable, dimension(:,:)   :: OBLdepth  !< Depth (positive) of OBL (m)
+  real, allocatable, dimension(:,:)   :: OBLdepth_original  !< Depth (positive) of OBL (m) without smoothing
   real, allocatable, dimension(:,:)   :: kOBL      !< Level (+fraction) of OBL extent
   real, allocatable, dimension(:,:)   :: OBLdepthprev !< previous Depth (positive) of OBL (m)
   real, allocatable, dimension(:,:,:) :: dRho      !< Bulk difference in density (kg/m3)
@@ -214,10 +217,16 @@ logical function KPP_init(paramFile, G, diag, Time, CS, passive, Waves)
                  'If False, calculates the non-local transport and tendencies but\n'//&
                  'purely for diagnostic purposes.',                                   &
                  default=.not. CS%passiveMode)
-  call get_param(paramFile, mdl, 'SMOOTH_BLD', CS%smoothBLD,  &
-                 'If True, applies a 1-1-4-1-1 Laplacian filter one time on HBLT.\n'//  &
-                 'computed via CVMix to reduce any horizontal two-grid-point noise.',   &
-                 default=.false.)
+  call get_param(paramFile, mdl, 'N_SMOOTH', CS%n_smooth,  &
+                 'The number of times the 1-1-4-1-1 Laplacian filter is applied on\n'//  &
+                 'OBL depth.',   &
+                 default=0)
+  if (CS%n_smooth > 0) then
+    call get_param(paramFile, mdl, 'DEEPEN_ONLY_VIA_SMOOTHING', CS%deepen_only,  &
+                   'If true, apply OBLdepth smoothing at a cell only if the OBLdepth.\n'// &
+                   'gets deeper via smoothing.',   &
+                   default=.false.)
+  endif
   call get_param(paramFile, mdl, 'RI_CRIT', CS%Ri_crit,                            &
                  'Critical bulk Richardson number used to define depth of the\n'// &
                  'surface Ocean Boundary Layer (OBL).',                            &
@@ -437,6 +446,10 @@ logical function KPP_init(paramFile, G, diag, Time, CS, passive, Waves)
       ! CMOR names are placeholders; must be modified by time period
       ! for CMOR compliance. Diag manager will be used for omlmax and
       ! omldamax.
+  CS%id_OBLdepth_original = register_diag_field('ocean_model', 'KPP_OBLdepth_original', diag%axesT1, Time, &
+      'Thickness of the surface Ocean Boundary Layer without smoothing calculated by [CVMix] KPP', 'meter', &
+      cmor_field_name='oml', cmor_long_name='ocean_mixed_layer_thickness_defined_by_mixing_scheme', &
+      cmor_units='m', cmor_standard_name='Ocean Mixed Layer Thickness Defined by Mixing Scheme')
   CS%id_BulkDrho = register_diag_field('ocean_model', 'KPP_BulkDrho', diag%axesTL, Time, &
       'Bulk difference in density used in Bulk Richardson number, as used by [CVMix] KPP', 'kg/m3')
   CS%id_BulkUz2 = register_diag_field('ocean_model', 'KPP_BulkUz2', diag%axesTL, Time, &
@@ -498,6 +511,7 @@ logical function KPP_init(paramFile, G, diag, Time, CS, passive, Waves)
   CS%OBLdepth(:,:) = 0.
   allocate( CS%kOBL( SZI_(G), SZJ_(G) ) )
   CS%kOBL(:,:) = 0.
+  if (CS%id_OBLdepth_original > 0) allocate( CS%OBLdepth_original( SZI_(G), SZJ_(G) ) )
 
   allocate( CS%OBLdepthprev( SZI_(G), SZJ_(G) ) );CS%OBLdepthprev(:,:)=0.0
   if (CS%id_BulkDrho > 0) allocate( CS%dRho( SZI_(G), SZJ_(G), SZK_(G) ) )
@@ -1063,6 +1077,7 @@ subroutine KPP_calculate(CS, G, GV, h, Temp, Salt, u, v, EOS, uStar, &
 
   ! send diagnostics to post_data
   if (CS%id_OBLdepth > 0) call post_data(CS%id_OBLdepth, CS%OBLdepth,        CS%diag)
+  if (CS%id_OBLdepth_original > 0) call post_data(CS%id_OBLdepth_original,CS%OBLdepth_original,CS%diag)
   if (CS%id_BulkDrho > 0) call post_data(CS%id_BulkDrho, CS%dRho,            CS%diag)
   if (CS%id_BulkUz2  > 0) call post_data(CS%id_BulkUz2,  CS%Uz2,             CS%diag)
   if (CS%id_BulkRi   > 0) call post_data(CS%id_BulkRi,   CS%BulkRi,          CS%diag)
@@ -1490,7 +1505,7 @@ subroutine KPP_compute_BLD(CS, G, GV, h, Temp, Salt, u, v, EOS, uStar, buoyFlux,
     enddo
   enddo
 
-  if (CS%smoothBLD) call KPP_smooth_BLD(CS,G,GV,h)
+  if (CS%n_smooth > 0) call KPP_smooth_BLD(CS,G,GV,h)
 
 end subroutine KPP_compute_BLD
 
@@ -1510,34 +1525,45 @@ subroutine KPP_smooth_BLD(CS,G,GV,h)
   real :: wc, ww, we, wn, ws ! averaging weights for smoothing
   real :: dh                 ! The local thickness used for calculating interface positions (m)
   real :: hcorr              ! A cumulative correction arising from inflation of vanished layers (m)
-  integer :: i, j, k
+  integer :: i, j, k, s
 
-  ! Update halos
-  call pass_var(CS%OBLdepth, G%Domain)
+  do s=1,CS%n_smooth
 
-  OBLdepth_original = CS%OBLdepth
+    ! Update halos
+    call pass_var(CS%OBLdepth, G%Domain)
 
-  ! apply smoothing on OBL depth
-  do j = G%jsc, G%jec
-    do i = G%isc, G%iec
+    OBLdepth_original = CS%OBLdepth
+    CS%OBLdepth_original = OBLdepth_original
 
-      ! skip land points
-      if (G%mask2dT(i,j)==0.) cycle
+    ! apply smoothing on OBL depth
+    do j = G%jsc, G%jec
+      do i = G%isc, G%iec
 
-      ! compute weights
-      ww = 0.125 * G%mask2dT(i-1,j)
-      we = 0.125 * G%mask2dT(i+1,j)
-      ws = 0.125 * G%mask2dT(i,j-1)
-      wn = 0.125 * G%mask2dT(i,j+1)
-      wc = 1.0 - (ww+we+wn+ws)
-
-      CS%OBLdepth(i,j) =  wc * OBLdepth_original(i,j)   &
-                        + ww * OBLdepth_original(i-1,j) &
-                        + we * OBLdepth_original(i+1,j) &
-                        + ws * OBLdepth_original(i,j-1) &
-                        + wn * OBLdepth_original(i,j+1)
+        ! skip land points
+        if (G%mask2dT(i,j)==0.) cycle
+
+        ! compute weights
+        ww = 0.125 * G%mask2dT(i-1,j)
+        we = 0.125 * G%mask2dT(i+1,j)
+        ws = 0.125 * G%mask2dT(i,j-1)
+        wn = 0.125 * G%mask2dT(i,j+1)
+        wc = 1.0 - (ww+we+wn+ws)
+
+        CS%OBLdepth(i,j) =  wc * OBLdepth_original(i,j)   &
+                          + ww * OBLdepth_original(i-1,j) &
+                          + we * OBLdepth_original(i+1,j) &
+                          + ws * OBLdepth_original(i,j-1) &
+                          + wn * OBLdepth_original(i,j+1)
+      enddo
     enddo
-  enddo
+
+    ! Apply OBLdepth smoothing at a cell only if the OBLdepth gets deeper via smoothing.
+    if (CS%deepen_only) CS%OBLdepth = max(CS%OBLdepth,CS%OBLdepth_original)
+
+    ! prevent OBL depths deeper than the bathymetric depth
+    where (CS%OBLdepth > G%bathyT) CS%OBLdepth = G%bathyT
+
+  enddo ! s-loop
 
   ! Update kOBL for smoothed OBL depths
   do j = G%jsc, G%jec
@@ -1559,7 +1585,7 @@ subroutine KPP_smooth_BLD(CS,G,GV,h)
         iFaceHeight(k+1) = iFaceHeight(k) - dh
       enddo
 
-      CS%kOBL(i,j)     = CVMix_kpp_compute_kOBL_depth( iFaceHeight, cellHeight, CS%OBLdepth(i,j) )
+      CS%kOBL(i,j) = CVMix_kpp_compute_kOBL_depth( iFaceHeight, cellHeight, CS%OBLdepth(i,j) )
 
     enddo
   enddo

src/parameterizations/vertical/MOM_set_diffusivity.F90

@@ -470,7 +470,7 @@ subroutine set_diffusivity(u, v, h, u_h, v_h, tv, fluxes, optics, visc, dt, &
 
     ! Add the Nikurashin and / or tidal bottom-driven mixing
     call calculate_tidal_mixing(h, N2_bot, j, TKE_to_Kd, maxTKE, G, GV, CS%tm_csp, &
-                                N2_lay, N2_int, Kd, Kd_int, CS%Kd_max)
+                                N2_lay, N2_int, Kd, Kd_int, CS%Kd_max, visc%Kv_slow)
 
     ! This adds the diffusion sustained by the energy extracted from the flow
     ! by the bottom drag.

src/parameterizations/vertical/MOM_tidal_mixing.F90

@@ -642,7 +642,7 @@ end function tidal_mixing_init
 !! tidal dissipation and to add the effect of internal-tide-driven mixing to the layer or interface
 !! diffusivities.
 subroutine calculate_tidal_mixing(h, N2_bot, j, TKE_to_Kd, max_TKE, G, GV, CS, &
-                                    N2_lay, N2_int, Kd, Kd_int, Kd_max)
+                                    N2_lay, N2_int, Kd, Kd_int, Kd_max, Kv)
   type(ocean_grid_type),                    intent(in)    :: G    !< The ocean's grid structure
   type(verticalGrid_type),                  intent(in)    :: GV   !< The ocean's vertical grid structure
   real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(in)    :: h    !< Layer thicknesses, in H (usually m or kg m-2)
@@ -655,10 +655,12 @@ subroutine calculate_tidal_mixing(h, N2_bot, j, TKE_to_Kd, max_TKE, G, GV, CS, &
   real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(inout) :: Kd
   real, dimension(SZI_(G),SZJ_(G),SZK_(G)+1), optional, intent(inout) :: Kd_int
   real,                                     intent(inout) :: Kd_max
+  real, dimension(:,:,:),                   pointer       :: Kv  !< The "slow" vertical viscosity at each interface
+                                                                 !! (not layer!) in m2 s-1.
 
   if (CS%Int_tide_dissipation .or. CS%Lee_wave_dissipation .or. CS%Lowmode_itidal_dissipation) then
     if (CS%use_CVMix_tidal) then
-      call calculate_CVMix_tidal(h, j, G, GV, CS, N2_int, Kd)
+      call calculate_CVMix_tidal(h, j, G, GV, CS, N2_int, Kd, Kv)
     else
       call add_int_tide_diffusivity(h, N2_bot, j, TKE_to_Kd, max_TKE, G, GV, CS, &
                                     N2_lay, Kd, Kd_int, Kd_max)
@@ -669,14 +671,16 @@ subroutine calculate_tidal_mixing(h, N2_bot, j, TKE_to_Kd, max_TKE, G, GV, CS, &
 
 !> Calls the CVMix routines to compute tidal dissipation and to add the effect of internal-tide-driven
 !! mixing to the interface diffusivities.
-subroutine calculate_CVMix_tidal(h, j, G, GV, CS, N2_int, Kd)
+subroutine calculate_CVMix_tidal(h, j, G, GV, CS, N2_int, Kd, Kv)
   integer,                                  intent(in)    :: j
   type(ocean_grid_type),                    intent(in)    :: G     !< Grid structure.
   type(verticalGrid_type),                  intent(in)    :: GV    !< ocean vertical grid structure
   type(tidal_mixing_cs),                    pointer       :: CS    !< This module's control structure.
   real, dimension(SZI_(G),SZK_(G)+1),       intent(in)    :: N2_int
   real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(in)    :: h     !< Layer thicknesses, in H (usually m or kg m-2).
   real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(inout) :: Kd
+  real, dimension(:,:,:),                   pointer       :: Kv    !< The "slow" vertical viscosity at each interface
+                                                                   !! (not layer!) in m2 s-1.
 
   ! local
   real, dimension(SZK_(G)+1) :: Kd_tidal    !< tidal diffusivity [m2/s]
@@ -744,11 +748,18 @@ subroutine calculate_CVMix_tidal(h, j, G, GV, CS, N2_int, Kd)
                                CVMix_params            = CS%CVMix_glb_params, &
                                CVMix_tidal_params_user = CS%CVMix_tidal_params)
 
+      ! Update diffusivity
       do k=1,G%ke
         Kd(i,j,k) = Kd(i,j,k) + 0.5*(Kd_tidal(k) + Kd_tidal(k+1) )
-        !TODO: Kv(i,j,k) = ????????????
       enddo
 
+      ! Update viscosity
+      if (associated(Kv)) then
+        do k=1,G%ke+1
+          Kv(i,j,k) = Kv(i,j,k) + Kv_tidal(k)
+        enddo
+      endif
+
       ! diagnostics
       if (associated(dd%Kd_itidal)) then
         dd%Kd_itidal(i,j,:) = Kd_tidal(:)
@@ -836,11 +847,18 @@ subroutine calculate_CVMix_tidal(h, j, G, GV, CS, N2_int, Kd)
                                           CVmix_params            = CS%CVMix_glb_params,  &
                                           CVmix_tidal_params_user = CS%CVMix_tidal_params)
 
+      ! Update diffusivity
       do k=1,G%ke
         Kd(i,j,k) = Kd(i,j,k) + 0.5*(Kd_tidal(k) + Kd_tidal(k+1) )
-        !TODO: Kv(i,j,k) = ????????????
       enddo
 
+      ! Update viscosity
+      if (associated(Kv)) then
+        do k=1,G%ke+1
+          Kv(i,j,k) = Kv(i,j,k) + Kv_tidal(k)
+        enddo
+      endif
+
       ! diagnostics
       if (associated(dd%Kd_itidal)) then
         dd%Kd_itidal(i,j,:) = Kd_tidal(:)
@@ -889,12 +907,7 @@ subroutine add_int_tide_diffusivity(h, N2_bot, j, TKE_to_Kd, max_TKE, G, GV, CS,
   real, dimension(SZI_(G),SZJ_(G),SZK_(G)+1), optional, intent(inout) :: Kd_int
   real,                                       intent(inout) :: Kd_max
 
-  ! This subroutine adds the effect of internal-tide-driven mixing to the layer diffusivities.
-  ! The mechanisms considered are (1) local dissipation of internal waves generated by the
-  ! barotropic flow ("itidal"), (2) local dissipation of internal waves generated by the propagating
-  ! low modes (rays) of the internal tide ("lowmode"), and (3) local dissipation of internal lee waves.
-  ! Will eventually need to add diffusivity due to other wave-breaking processes (e.g. Bottom friction,
-  ! Froude-number-depending breaking, PSI, etc.).
+  ! local
 
   real, dimension(SZI_(G)) :: &
     htot,             & ! total thickness above or below a layer, or the