Merge pull request #106 from gustavo-marques/bobs_solution

Prohibits USE_GME to be used when SPLIT=False

.gitlab-ci.yml

@@ -32,6 +32,7 @@ setup:
     - git clone --recursive http://gitlab.gfdl.noaa.gov/ogrp/Gaea-stats-MOM6-examples.git tests && cd tests
     # Install / update testing scripts
     - git clone https://github.com/adcroft/MRS.git MRS
+    - (cd MRS ; git checkout xanadu-fms)
     # Update MOM6-examples and submodules
     - (cd MOM6-examples && git checkout . && git checkout dev/gfdl && git pull && git submodule init && git submodule update)
     - (cd MOM6-examples/src/MOM6 && git submodule update)
@@ -73,7 +74,7 @@ gnu:ice-ocean-nolibs:
     - time tar zxf $CACHE_DIR/tests_$CI_PIPELINE_ID.tgz && cd tests
     - make -f MRS/Makefile.build build/gnu/env && cd build/gnu
     # mkdir -p build/gnu/repro/symmetric_dynamic/ocean_only && cd build/gnu/repro/symmetric_dynamic/ocean_only
-    - ../../MOM6-examples/src/mkmf/bin/list_paths -l ../../../config_src/{coupled_driver,dynamic} ../../../src ../../MOM6-examples/src/{FMS,coupler,SIS2,icebergs,ice_ocean_extras,land_null,atmos_null}
+    - ../../MOM6-examples/src/mkmf/bin/list_paths -l ../../../config_src/{coupled_driver,dynamic} ../../../src ../../MOM6-examples/src/{FMS,coupler,SIS2,icebergs,ice_param,land_null,atmos_null}
     - ../../MOM6-examples/src/mkmf/bin/mkmf -t ../../MOM6-examples/src/mkmf/templates/ncrc-gnu.mk -p MOM6 -c"-Duse_libMPI -Duse_netCDF -D_USE_LEGACY_LAND_ -Duse_AM3_physics" path_names
     - time (source ./env ; make NETCDF=3 REPRO=1 MOM6 -s -j)
 

.testing/Makefile

@@ -18,7 +18,7 @@ else
 EXPERIMENTS ?= unit_tests double_gyre flow_downslope/z CVmix_SCM_tests/cooling_only/EPBL
 endif
 
-FMS_PACKAGES ?= platform,include,memutils,constants,mpp,fms,time_manager,diag_manager,data_override,coupler/ensemble_manager.F90,axis_utils,horiz_interp,time_interp,astronomy,mosaic,random_numbers
+FMS_PACKAGES ?= platform,include,memutils,constants,mpp,fms,time_manager,diag_manager,data_override,coupler/coupler_types.F90,coupler/ensemble_manager.F90,axis_utils,horiz_interp,time_interp,astronomy,mosaic,random_numbers
 TEMPLATE ?= .testing/linux-ubuntu-xenial-gnu.mk
 MPIRUN ?= mpirun
 

src/core/MOM.F90

@@ -999,13 +999,11 @@ subroutine step_MOM_dynamics(forces, p_surf_begin, p_surf_end, dt, dt_thermo, &
     if (CS%use_RK2) then
       call step_MOM_dyn_unsplit_RK2(u, v, h, CS%tv, CS%visc, Time_local, dt, forces, &
                p_surf_begin, p_surf_end, CS%uh, CS%vh, CS%uhtr, CS%vhtr, &
-               CS%eta_av_bc, G, GV, US, CS%dyn_unsplit_RK2_CSp, CS%VarMix, CS%MEKE, &
-               CS%Barotropic_CSp, CS%thickness_diffuse_CSp)
+               CS%eta_av_bc, G, GV, US, CS%dyn_unsplit_RK2_CSp, CS%VarMix, CS%MEKE)
     else
       call step_MOM_dyn_unsplit(u, v, h, CS%tv, CS%visc, Time_local, dt, forces, &
                p_surf_begin, p_surf_end, CS%uh, CS%vh, CS%uhtr, CS%vhtr, &
-               CS%eta_av_bc, G, GV, US, CS%dyn_unsplit_CSp, CS%VarMix, CS%MEKE, &
-               CS%Barotropic_CSp, CS%thickness_diffuse_CSp,  Waves=Waves)
+               CS%eta_av_bc, G, GV, US, CS%dyn_unsplit_CSp, CS%VarMix, CS%MEKE, Waves=Waves)
     endif
     if (showCallTree) call callTree_waypoint("finished step_MOM_dyn_unsplit (step_MOM)")
 

src/core/MOM_dynamics_split_RK2.F90

@@ -689,8 +689,8 @@ subroutine step_MOM_dyn_split_RK2(u, v, h, tv, visc, &
 ! diffu = horizontal viscosity terms (u_av)
   call cpu_clock_begin(id_clock_horvisc)
   call horizontal_viscosity(u_av, v_av, h_av, CS%diffu, CS%diffv, &
-                            MEKE, Varmix, CS%barotropic_CSp, thickness_diffuse_CSp, &
-                            G, GV, US, CS%hor_visc_CSp, OBC=CS%OBC)
+                            MEKE, Varmix, G, GV, US, CS%hor_visc_CSp, &
+                            OBC=CS%OBC, BT=CS%barotropic_CSp)
   call cpu_clock_end(id_clock_horvisc)
   if (showCallTree) call callTree_wayPoint("done with horizontal_viscosity (step_MOM_dyn_split_RK2)")
 
@@ -1147,14 +1147,11 @@ subroutine initialize_dyn_split_RK2(u, v, h, uh, vh, eta, Time, G, GV, US, param
                        CS%barotropic_CSp, restart_CS, calc_dtbt, CS%BT_cont, &
                        CS%tides_CSp)
 
-!  CS%barotropic_CSp => Barotropic_CSp
-!  CS%thickness_diffuse_CSp => thickness_diffuse_CSp
-
   if (.not. query_initialized(CS%diffu,"diffu",restart_CS) .or. &
       .not. query_initialized(CS%diffv,"diffv",restart_CS)) &
     call horizontal_viscosity(u, v, h, CS%diffu, CS%diffv, MEKE, VarMix, &
-                              CS%barotropic_CSp, thickness_diffuse_CSp, &
-                              G, GV, US, CS%hor_visc_CSp, OBC=CS%OBC)
+                              G, GV, US, CS%hor_visc_CSp, &
+                              OBC=CS%OBC, BT=CS%barotropic_CSp)
   if (.not. query_initialized(CS%u_av,"u2", restart_CS) .or. &
       .not. query_initialized(CS%u_av,"v2", restart_CS)) then
     CS%u_av(:,:,:) = u(:,:,:)

src/core/MOM_dynamics_unsplit.F90

@@ -183,7 +183,7 @@ module MOM_dynamics_unsplit
 !! 3rd order (for the inviscid momentum equations) order scheme
 subroutine step_MOM_dyn_unsplit(u, v, h, tv, visc, Time_local, dt, forces, &
                   p_surf_begin, p_surf_end, uh, vh, uhtr, vhtr, eta_av, G, GV, US, CS, &
-                  VarMix, MEKE, Barotropic, thickness_diffuse, Waves)
+                  VarMix, MEKE, Waves)
   type(ocean_grid_type),   intent(inout) :: 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
@@ -218,10 +218,6 @@ subroutine step_MOM_dyn_unsplit(u, v, h, tv, visc, Time_local, dt, forces, &
                                  !! that specify the spatially variable viscosities.
   type(MEKE_type),         pointer       :: MEKE   !< A pointer to a structure containing
                                  !! fields related to the Mesoscale Eddy Kinetic Energy.
-  type(barotropic_CS),      pointer       :: Barotropic!< Pointer to a structure containing
-                                                              !! barotropic velocities
-  type(thickness_diffuse_CS),  pointer       :: thickness_diffuse!< Pointer to a structure containing
-                                                              !! interface height diffusivities
   type(wave_parameters_CS), optional, pointer :: Waves !< A pointer to a structure containing
                                  !! fields related to the surface wave conditions
 
@@ -260,7 +256,7 @@ subroutine step_MOM_dyn_unsplit(u, v, h, tv, visc, Time_local, dt, forces, &
   call enable_averaging(dt,Time_local, CS%diag)
   call cpu_clock_begin(id_clock_horvisc)
   call horizontal_viscosity(u, v, h, CS%diffu, CS%diffv, MEKE, Varmix, &
-                            Barotropic, thickness_diffuse, G, GV, US, CS%hor_visc_CSp)
+                            G, GV, US, CS%hor_visc_CSp)
   call cpu_clock_end(id_clock_horvisc)
   call disable_averaging(CS%diag)
 

src/core/MOM_dynamics_unsplit_RK2.F90

@@ -186,7 +186,7 @@ module MOM_dynamics_unsplit_RK2
 !> Step the MOM6 dynamics using an unsplit quasi-2nd order Runge-Kutta scheme
 subroutine step_MOM_dyn_unsplit_RK2(u_in, v_in, h_in, tv, visc, Time_local, dt, forces, &
                   p_surf_begin, p_surf_end, uh, vh, uhtr, vhtr, eta_av, G, GV, US, CS, &
-                  VarMix, MEKE, Barotropic, thickness_diffuse)
+                  VarMix, MEKE)
   type(ocean_grid_type),             intent(inout) :: 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
@@ -232,10 +232,6 @@ subroutine step_MOM_dyn_unsplit_RK2(u_in, v_in, h_in, tv, visc, Time_local, dt,
   type(MEKE_type),                   pointer       :: MEKE    !< A pointer to a structure containing
                                                               !! fields related to the Mesoscale
                                                               !! Eddy Kinetic Energy.
-  type(barotropic_CS),               pointer       :: Barotropic!< Pointer to a structure containing
-                                                              !! barotropic velocities
-  type(thickness_diffuse_CS),        pointer       :: thickness_diffuse!< Pointer to a structure containing
-                                                              !! interface height diffusivities
   ! Local variables
   real, dimension(SZI_(G),SZJ_(G),SZK_(G)) :: h_av, hp
   real, dimension(SZIB_(G),SZJ_(G),SZK_(G)) :: up
@@ -271,7 +267,7 @@ subroutine step_MOM_dyn_unsplit_RK2(u_in, v_in, h_in, tv, visc, Time_local, dt,
   call enable_averaging(dt,Time_local, CS%diag)
   call cpu_clock_begin(id_clock_horvisc)
   call horizontal_viscosity(u_in, v_in, h_in, CS%diffu, CS%diffv, MEKE, VarMix, &
-                            Barotropic, thickness_diffuse, G, GV, US, CS%hor_visc_CSp)
+                            G, GV, US, CS%hor_visc_CSp)
   call cpu_clock_end(id_clock_horvisc)
   call disable_averaging(CS%diag)
   call pass_vector(CS%diffu, CS%diffv, G%Domain, clock=id_clock_pass)

src/diagnostics/MOM_obsolete_params.F90

@@ -69,6 +69,9 @@ subroutine find_obsolete_params(param_file)
        hint="Instead use OBC_SEGMENT_XXX_DATA.")
   call obsolete_char(param_file, "EXTEND_OBC_SEGMENTS", &
        hint="This option is no longer needed, nor supported.")
+  call obsolete_char(param_file, "MEKE_VISCOSITY_COEFF", &
+       hint="This option has been replaced by MEKE_VISCOSITY_COEFF_KU and \n" //&
+            " MEKE_VISCOSITY_COEFF_AU. Please set these parameters instead.")
   nseg = 0
   call read_param(param_file, "OBC_NUMBER_OF_SEGMENTS", nseg)
   do l_seg = 1,nseg

src/parameterizations/lateral/MOM_MEKE.F90

@@ -1420,9 +1420,13 @@ end subroutine MEKE_end
 !! \subsection section_MEKE_viscosity Viscosity derived from MEKE
 !!
 !! As for \f$ \kappa_M \f$, the predicted eddy velocity scale can be
-!! used to form an eddy viscosity:
+!! used to form a harmonic eddy viscosity,
 !!
-!! \f[  \kappa_u = \gamma_u \sqrt{ U_e^2 A_\Delta } . \f]
+!! \f[  \kappa_u = \gamma_u \sqrt{ U_e^2 A_\Delta }  \f]
+!!
+!! as well as a biharmonic eddy viscosity,
+!!
+!! \f[  \kappa_4 = \gamma_4 \sqrt{ U_e^2 A_\Delta^3 }  \f]
 !!
 !! \subsection section_MEKE_limit_case Limit cases for local source-dissipative balance
 !!
@@ -1459,7 +1463,8 @@ end subroutine MEKE_end
 !! | \f$ \kappa_4 \f$      | <code>MEKE_K4</code> |
 !! | \f$ \gamma_\kappa \f$ | <code>MEKE_KHCOEFF</code> |
 !! | \f$ \gamma_M \f$      | <code>MEKE_KHMEKE_FAC</code> |
-!! | \f$ \gamma_u \f$      | <code>MEKE_VISCOSITY_COEFF</code> |
+!! | \f$ \gamma_u \f$      | <code>MEKE_VISCOSITY_COEFF_KU</code> |
+!! | \f$ \gamma_4 \f$      | <code>MEKE_VISCOSITY_COEFF_AU</code> |
 !! | \f$ \gamma_{min}^2 \f$| <code>MEKE_MIN_GAMMA2</code> |
 !! | \f$ \alpha_d \f$      | <code>MEKE_ALPHA_DEFORM</code> |
 !! | \f$ \alpha_f \f$      | <code>MEKE_ALPHA_FRICT</code> |

src/parameterizations/lateral/MOM_hor_visc.F90

@@ -190,16 +190,16 @@ module MOM_hor_visc
 !!   u[is-2:ie+2,js-2:je+2]
 !!   v[is-2:ie+2,js-2:je+2]
 !!   h[is-1:ie+1,js-1:je+1]
-subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, Barotropic, &
-                                thickness_diffuse, G, GV, US, CS, OBC)
+subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US, &
+                                CS, OBC, BT)
   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(SZIB_(G),SZJ_(G),SZK_(G)), &
                                  intent(in)  :: u      !< The zonal velocity [m s-1].
   real, dimension(SZI_(G),SZJB_(G),SZK_(G)), &
                                  intent(in)  :: v      !< The meridional velocity [m s-1].
   real, dimension(SZI_(G),SZJ_(G),SZK_(G)),  &
-                                 intent(inout)  :: h      !< Layer thicknesses, in H
+                                 intent(inout)  :: h   !< Layer thicknesses, in H
                                                        !! (usually m or kg m-2).
   real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), &
                                  intent(out) :: diffu  !< Zonal acceleration due to convergence of
@@ -211,14 +211,12 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, Barotropic,
                                                        !! related to Mesoscale Eddy Kinetic Energy.
   type(VarMix_CS),               pointer     :: VarMix !< Pointer to a structure with fields that
                                                        !! specify the spatially variable viscosities
-  type(barotropic_CS),           pointer     :: Barotropic  !< Pointer to a structure containing
-                                                       !! barotropic velocities
-  type(thickness_diffuse_CS),    pointer     :: thickness_diffuse  !< Pointer to a structure containing
-                                                       !! interface height diffusivities
   type(hor_visc_CS),             pointer     :: CS     !< Control structure returned by a previous
   type(unit_scale_type),         intent(in)  :: US     !< A dimensional unit scaling type
                                                        !! call to hor_visc_init.
   type(ocean_OBC_type), optional, pointer    :: OBC    !< Pointer to an open boundary condition type
+  type(barotropic_CS),  optional, pointer    :: BT     !< Pointer to a structure containing
+                                                       !! barotropic velocities.
 
   ! Local variables
   real, dimension(SZIB_(G),SZJ_(G)) :: &
@@ -386,20 +384,6 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, Barotropic,
   use_MEKE_Ku = associated(MEKE%Ku)
   use_MEKE_Au = associated(MEKE%Au)
 
-!$OMP parallel do default(none) shared(Isq,Ieq,Jsq,Jeq,nz,CS,G,GV,u,v,is,js,ie,je,h,  &
-!$OMP                                  rescale_Kh,VarMix,h_neglect,h_neglect3,        &
-!$OMP                                  Kh_h,Ah_h,Kh_q,Ah_q,diffu,apply_OBC,OBC,diffv, &
-!$OMP                                  find_FrictWork,FrictWork,use_MEKE_Ku,          &
-!$OMP                                  use_MEKE_Au, MEKE, hq,                         &
-!$OMP                                  mod_Leith, legacy_bound, div_xx_h, vort_xy_q)  &
-!$OMP                          private(u0, v0, sh_xx, str_xx, visc_bound_rem, &
-!$OMP                                  sh_xy, str_xy, Ah, Kh, AhSm, KhSm, dvdx, dudy, &
-!$OMP                                  sh_xx_bt, sh_xy_bt, dvdx_bt, dudy_bt, &
-!$OMP                                  bhstr_xx, bhstr_xy,FatH,RoScl, hu, hv, h_u, h_v, &
-!$OMP                                  vort_xy,vort_xy_dx,vort_xy_dy,Vort_mag,AhLth,KhLth, &
-!$OMP                                  div_xx, div_xx_dx, div_xx_dy,local_strain,          &
-!$OMP                                  Shear_mag, h2uq, h2vq, Kh_scale, hrat_min)
-
   do j=js,je ; do i=is,ie
     boundary_mask(i,j) = (G%mask2dCu(I,j) * G%mask2dCv(i,J) * G%mask2dCu(I-1,j) * G%mask2dCv(i,J-1))
   enddo ; enddo
@@ -419,7 +403,7 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, Barotropic,
 !    call pass_var(boundary_mask, G%Domain, complete=.true.)
 
     ! Get barotropic velocities and their gradients
-    call barotropic_get_tav(Barotropic, ubtav, vbtav, G)
+    call barotropic_get_tav(BT, ubtav, vbtav, G)
     call pass_vector(ubtav, vbtav, G%Domain)
 
     do j=js,je ; do i=is,ie
@@ -480,6 +464,19 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, Barotropic,
 
   endif ! use_GME
 
+  !$OMP parallel do default(none) shared(Isq,Ieq,Jsq,Jeq,nz,CS,G,GV,u,v,is,js,ie,je,h,  &
+  !$OMP                                  rescale_Kh,VarMix,h_neglect,h_neglect3,        &
+  !$OMP                                  Kh_h,Ah_h,Kh_q,Ah_q,diffu,apply_OBC,OBC,diffv, &
+  !$OMP                                  find_FrictWork,FrictWork,use_MEKE_Ku,          &
+  !$OMP                                  use_MEKE_Au, MEKE, hq,                         &
+  !$OMP                                  mod_Leith, legacy_bound, div_xx_h, vort_xy_q)  &
+  !$OMP                          private(u0, v0, sh_xx, str_xx, visc_bound_rem, &
+  !$OMP                                  sh_xy, str_xy, Ah, Kh, AhSm, KhSm, dvdx, dudy, &
+  !$OMP                                  sh_xx_bt, sh_xy_bt, dvdx_bt, dudy_bt, &
+  !$OMP                                  bhstr_xx, bhstr_xy,FatH,RoScl, hu, hv, h_u, h_v, &
+  !$OMP                                  vort_xy,vort_xy_dx,vort_xy_dy,Vort_mag,AhLth,KhLth, &
+  !$OMP                                  div_xx, div_xx_dx, div_xx_dy,local_strain,          &
+  !$OMP                                  Shear_mag, h2uq, h2vq, Kh_scale, hrat_min)
   do k=1,nz
 
     ! The following are the forms of the horizontal tension and horizontal
@@ -1410,6 +1407,8 @@ subroutine hor_visc_init(Time, G, US, param_file, diag, CS)
   logical :: get_all       ! If true, read and log all parameters, regardless of
                            ! whether they are used, to enable spell-checking of
                            ! valid parameters.
+  logical :: split         ! If true, use the split time stepping scheme.
+                           ! If false and USE_GME = True, issue a FATAL error.
   character(len=64) :: inputdir, filename
   real    :: deg2rad       ! Converts degrees to radians
   real    :: slat_fn       ! sin(lat)**Kh_pwr_of_sine
@@ -1563,7 +1562,7 @@ subroutine hor_visc_init(Time, G, US, param_file, diag, CS)
                  "the grid spacing to calculate the biharmonic viscosity. \n"//&
                  "The final viscosity is the largest of this scaled \n"//&
                  "viscosity, the Smagorinsky and Leith viscosities, and AH.", &
-                 units="m s-1", default=0.1)
+                 units="m s-1", default=0.0)
     call get_param(param_file, mdl, "SMAGORINSKY_AH", CS%Smagorinsky_Ah, &
                  "If true, use a biharmonic Smagorinsky nonlinear eddy \n"//&
                  "viscosity.", default=.false.)
@@ -1642,6 +1641,14 @@ subroutine hor_visc_init(Time, G, US, param_file, diag, CS)
                  "If true, use the GM+E backscatter scheme in association \n"//&
                  "with the Gent and McWilliams parameterization.", default=.false.)
 
+  if (CS%use_GME) then
+    call get_param(param_file, mdl, "SPLIT", split, &
+                 "Use the split time stepping if true.", default=.true., &
+                  do_not_log=.true.)
+    if (.not. split) call MOM_error(FATAL,"ERROR: Currently, USE_GME = True "// &
+                                           "cannot be used with SPLIT=False.")
+  endif
+
   if (CS%bound_Kh .or. CS%bound_Ah .or. CS%better_bound_Kh .or. CS%better_bound_Ah) &
     call get_param(param_file, mdl, "DT", dt, &
                  "The (baroclinic) dynamics time step.", units = "s", &