Merge branch 'dev/master' of github.com:NOAA-GFDL/MOM6 into dev/master

src/core/MOM.F90

@@ -1424,36 +1424,6 @@ subroutine initialize_MOM(Time, param_file, dirs, CS, Time_in)
 
   call find_obsolete_params(param_file)
 
-#ifdef SYMMETRIC_MEMORY_
-  symmetric = .true.
-#else
-  symmetric = .false.
-#endif
-#ifdef STATIC_MEMORY_
-  call MOM_domains_init(G%domain, param_file, symmetric=symmetric, &
-            static_memory=.true., NIHALO=NIHALO_, NJHALO=NJHALO_, &
-            NIGLOBAL=NIGLOBAL_, NJGLOBAL=NJGLOBAL_, NIPROC=NIPROC_, &
-            NJPROC=NJPROC_)
-#else
-  call MOM_domains_init(G%domain, param_file, symmetric=symmetric)
-#endif
-  call callTree_waypoint("domains initialized (initialize_MOM)")
-
-  call MOM_checksums_init(param_file)
-
-  call diag_mediator_infrastructure_init()
-  call MOM_io_init(param_file)
-  call MOM_grid_init(G, param_file)
-
-  call create_dyn_horgrid(dG, G%HI)
-  dG%first_direction = G%first_direction
-  dG%bathymetry_at_vel = G%bathymetry_at_vel
-  call clone_MOM_domain(G%Domain, dG%Domain)
-
-  call verticalGridInit( param_file, CS%GV )
-  GV => CS%GV
-  dG%g_Earth = GV%g_Earth
-
   ! Read relevant parameters and write them to the model log.
   call log_version(param_file, "MOM", version, "")
   call get_param(param_file, "MOM", "VERBOSITY", verbosity,  &
@@ -1663,6 +1633,38 @@ subroutine initialize_MOM(Time, param_file, dirs, CS, Time_in)
   if (CS%adiabatic .and. CS%bulkmixedlayer) call MOM_error(FATAL, &
     "MOM: ADIABATIC and BULKMIXEDLAYER can not both be defined.")
 
+  ! Set up the model domain and grids.
+#ifdef SYMMETRIC_MEMORY_
+  symmetric = .true.
+#else
+  symmetric = .false.
+#endif
+#ifdef STATIC_MEMORY_
+  call MOM_domains_init(G%domain, param_file, symmetric=symmetric, &
+            static_memory=.true., NIHALO=NIHALO_, NJHALO=NJHALO_, &
+            NIGLOBAL=NIGLOBAL_, NJGLOBAL=NJGLOBAL_, NIPROC=NIPROC_, &
+            NJPROC=NJPROC_)
+#else
+  call MOM_domains_init(G%domain, param_file, symmetric=symmetric)
+#endif
+  call callTree_waypoint("domains initialized (initialize_MOM)")
+
+  call MOM_checksums_init(param_file)
+
+  call diag_mediator_infrastructure_init()
+  call MOM_io_init(param_file)
+  call MOM_grid_init(G, param_file)
+
+  call create_dyn_horgrid(dG, G%HI)
+  dG%first_direction = G%first_direction
+  dG%bathymetry_at_vel = G%bathymetry_at_vel
+  call clone_MOM_domain(G%Domain, dG%Domain)
+
+  call verticalGridInit( param_file, CS%GV )
+  GV => CS%GV
+  dG%g_Earth = GV%g_Earth
+
+
   ! Allocate the auxiliary non-symmetric domain for debugging or I/O purposes.
   if (CS%debug .or. dG%symmetric) &
     call clone_MOM_domain(dG%Domain, dG%Domain_aux, symmetric=.false.)

src/core/MOM_CoriolisAdv.F90

@@ -875,7 +875,7 @@ subroutine CoriolisAdv_init(Time, G, param_file, diag, AD, CS)
   CS%diag => diag ; CS%Time => Time
 
   ! Read all relevant parameters and write them to the model log.
-  call log_version(param_file, mod, version)
+  call log_version(param_file, mod, version, "")
   call get_param(param_file, mod, "NOSLIP", CS%no_slip, &
                  "If true, no slip boundary conditions are used; otherwise \n"//&
                  "free slip boundary conditions are assumed. The \n"//&

src/core/MOM_PressureForce.F90

@@ -144,7 +144,7 @@ subroutine PressureForce_init(Time, G, GV, param_file, diag, CS, tides_CSp)
   else ; allocate(CS) ; endif
 
   ! Read all relevant parameters and write them to the model log.
-  call log_version(param_file, mod, version)
+  call log_version(param_file, mod, version, "")
   call get_param(param_file, mod, "ANALYTIC_FV_PGF", CS%Analytic_FV_PGF, &
                  "If true the pressure gradient forces are calculated \n"//&
                  "with a finite volume form that analytically integrates \n"//&

src/core/MOM_PressureForce_Montgomery.F90

@@ -967,7 +967,7 @@ subroutine PressureForce_Mont_init(Time, G, GV, param_file, diag, CS, tides_CSp)
   endif
 
   mod = "MOM_PressureForce_Mont"
-  call log_version(param_file, mod, version)
+  call log_version(param_file, mod, version, "")
   call get_param(param_file, mod, "RHO_0", CS%Rho0, &
                  "The mean ocean density used with BOUSSINESQ true to \n"//&
                  "calculate accelerations and the mass for conservation \n"//&

src/core/MOM_PressureForce_analytic_FV.F90

@@ -871,7 +871,7 @@ subroutine PressureForce_AFV_init(Time, G, GV, param_file, diag, CS, tides_CSp)
   endif
 
   mod = "MOM_PressureForce_AFV"
-  call log_version(param_file, mod, version)
+  call log_version(param_file, mod, version, "")
   call get_param(param_file, mod, "RHO_0", CS%Rho0, &
                  "The mean ocean density used with BOUSSINESQ true to \n"//&
                  "calculate accelerations and the mass for conservation \n"//&

src/core/MOM_continuity.F90

@@ -190,7 +190,7 @@ subroutine continuity_init(Time, G, GV, param_file, diag, CS)
   allocate(CS)
 
   ! Read all relevant parameters and write them to the model log.
-  call log_version(param_file, mod, version)
+  call log_version(param_file, mod, version, "")
   call get_param(param_file, mod, "CONTINUITY_SCHEME", tmpstr, &
                  "CONTINUITY_SCHEME selects the discretization for the \n"//&
                  "continuity solver. The only valid value currently is: \n"//&

src/core/MOM_continuity_PPM.F90

@@ -2135,7 +2135,7 @@ subroutine continuity_PPM_init(Time, G, GV, param_file, diag, CS)
   allocate(CS)
 
 ! Read all relevant parameters and write them to the model log.
-  call log_version(param_file, mod, version)
+  call log_version(param_file, mod, version, "")
   call get_param(param_file, mod, "MONOTONIC_CONTINUITY", CS%monotonic, &
                  "If true, CONTINUITY_PPM uses the Colella and Woodward \n"//&
                  "monotonic limiter.  The default (false) is to use a \n"//&

src/core/MOM_open_boundary.F90

@@ -220,7 +220,7 @@ subroutine open_boundary_init(Time, G, param_file, diag, CS)
     return
   endif
 
-  call log_version(param_file, mod, version)
+  call log_version(param_file, mod, version, "")
   call get_param(param_file, mod, "APPLY_OBC_U_FLATHER_EAST", flather_east, &
                  "If true, some zonal velocity points use Flather open \n"//&
                  "boundary conditions on the east side of the ocean.", &

src/diagnostics/MOM_diag_to_Z.F90

@@ -988,7 +988,7 @@ subroutine MOM_diag_to_Z_init(Time, G, GV, param_file, diag, CS)
   CS%diag => diag
 
   ! Read parameters and write them to the model log.
-  call log_version(param_file, mod, version)
+  call log_version(param_file, mod, version, "")
   ! Read in z-space info from a NetCDF file.
   call get_param(param_file, mod, "Z_OUTPUT_GRID_FILE", zgrid_file, &
                  "The file that specifies the vertical grid for \n"//&

src/equation_of_state/MOM_EOS.F90

@@ -446,7 +446,7 @@ subroutine EOS_init(param_file, EOS)
   if (.not.associated(EOS)) call EOS_allocate(EOS)
 
   ! Read all relevant parameters and write them to the model log.
-  call log_version(param_file, mod, version)
+  call log_version(param_file, mod, version, "")
 
   call get_param(param_file, mod, "EQN_OF_STATE", tmpstr, &
                  "EQN_OF_STATE determines which ocean equation of state \n"//&

src/framework/MOM_diag_mediator.F90

@@ -215,7 +215,7 @@ subroutine set_axes_info(G, GV, param_file, diag_cs, set_vertical)
   set_vert = .true. ; if (present(set_vertical)) set_vert = set_vertical
 
   ! Read all relevant parameters and write them to the model log.
-  call log_version(param_file, mod, version)
+  call log_version(param_file, mod, version, "")
 
   if(G%symmetric) then
     id_xq = diag_axis_init('xq', G%gridLonB(G%isgB:G%iegB), G%x_axis_units, 'x', &

src/framework/MOM_domains.F90

@@ -941,7 +941,7 @@ subroutine MOM_domains_init(MOM_dom, param_file, symmetric, static_memory, &
   endif
 
   ! Read all relevant parameters and write them to the model log.
-  call log_version(param_file, mod, version)
+  call log_version(param_file, mod, version, "")
   call get_param(param_file, mod, "REENTRANT_X", reentrant_x, &
                  "If true, the domain is zonally reentrant.", default=.true.)
   call get_param(param_file, mod, "REENTRANT_Y", reentrant_y, &

src/framework/MOM_file_parser.F90

@@ -284,7 +284,7 @@ subroutine close_param_file(CS, quiet_close, component)
 
   ! Log the parameters for the parser.
   mod = "MOM_file_parser"
-  call log_version(CS, mod, version)
+  call log_version(CS, mod, version, "")
   call log_param(CS, mod, "SEND_LOG_TO_STDOUT", &
                         CS%log_to_stdout, &
                  "If true, all log messages are also sent to stdout.", &

src/framework/MOM_restart.F90

@@ -1237,7 +1237,7 @@ subroutine restart_init(param_file, CS, restart_root)
   allocate(CS)
 
   ! Read all relevant parameters and write them to the model log.
-  call log_version(param_file, mod, version)
+  call log_version(param_file, mod, version, "")
   call get_param(param_file, mod, "PARALLEL_RESTARTFILES", &
                                 CS%parallel_restartfiles, &
                  "If true, each processor writes its own restart file, \n"//&

src/framework/MOM_write_cputime.F90

@@ -100,7 +100,7 @@ subroutine MOM_write_cputime_init(param_file, directory, Input_start_time, CS)
   endif
 
   ! Read all relevant parameters and write them to the model log.
-  call log_version(param_file, mod, version)
+  call log_version(param_file, mod, version, "")
   call get_param(param_file, mod, "MAXCPU", CS%maxcpu, &
                  "The maximum amount of cpu time per processor for which \n"//&
                  "MOM should run before saving a restart file and \n"//&

src/initialization/MOM_coord_initialization.F90

@@ -53,7 +53,7 @@ subroutine MOM_initialize_coord(GV, PF, write_geom, output_dir, tv, max_depth)
   nz = GV%ke
 
   call callTree_enter("MOM_initialize_coord(), MOM_coord_initialization.F90")
-!  call log_version(PF, mod, version)
+  call log_version(PF, mod, version, "")
   call get_param(PF, mod, "DEBUG", debug, default=.false.)
 
 ! Set-up the layer densities, GV%Rlay, and reduced gravities, GV%g_prime.

src/initialization/MOM_fixed_initialization.F90

@@ -59,7 +59,7 @@ subroutine MOM_initialize_fixed(G, PF, write_geom, output_dir)
 #include "version_variable.h"
 
   call callTree_enter("MOM_initialize_fixed(), MOM_fixed_initialization.F90")
-  call log_version(PF, mod, version)
+  call log_version(PF, mod, version, "")
   call get_param(PF, mod, "DEBUG", debug, default=.false.)
 
   call get_param(PF, mod, "INPUTDIR", inputdir, &

src/initialization/MOM_state_initialization.F90

@@ -159,7 +159,7 @@ subroutine MOM_initialize_state(u, v, h, tv, Time, G, GV, PF, dirs, &
   IsdB = G%IsdB ; IedB = G%IedB ; JsdB = G%JsdB ; JedB = G%JedB
 
   call callTree_enter("MOM_initialize_state(), MOM_state_initialization.F90")
-  call log_version(PF, mod, version)
+  call log_version(PF, mod, version, "")
   call get_param(PF, mod, "DEBUG", debug, default=.false.)
 
   new_sim = .false.
@@ -2123,7 +2123,6 @@ subroutine MOM_temp_salt_initialize_from_Z(h, tv, G, GV, PF, dirs)
 !                      model parameter values.
 !  (in)      dirs    - A structure containing several relevant directory paths.
 
-
   type(ocean_grid_type),                 intent(inout) :: G
   real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(out)   :: h    
   type(thermo_var_ptrs),                 intent(inout) :: tv
@@ -2140,12 +2139,8 @@ subroutine MOM_temp_salt_initialize_from_Z(h, tv, G, GV, PF, dirs)
 
 ! This include declares and sets the variable "version".
 #include "version_variable.h"
-
-
   character(len=40)  :: mod = "MOM_initialize_layers_from_Z" ! This module's name.
 
-
-
   integer :: is, ie, js, je, nz ! compute domain indices
   integer :: isc,iec,jsc,jec    ! global compute domain indices
   integer :: isg, ieg, jsg, jeg ! global extent
@@ -2208,7 +2203,7 @@ subroutine MOM_temp_salt_initialize_from_Z(h, tv, G, GV, PF, dirs)
   PI_180=atan(1.0)/45.
 
   call callTree_enter(trim(mod)//"(), MOM_state_initialization.F90")
-  call log_version(PF, mod, version)
+  call log_version(PF, mod, version, "")
 
   new_sim = .false.
   if ((dirs%input_filename(1:1) == 'n') .and. &

src/parameterizations/vertical/MOM_ALE_sponge.F90

@@ -93,7 +93,7 @@ subroutine initialize_ALE_sponge(Iresttime, data_h, nz_data, G, param_file, CS)
   endif
 
 ! Set default, read and log parameters
-  call log_version(param_file, mod, version)
+  call log_version(param_file, mod, version, "")
   call get_param(param_file, mod, "SPONGE", use_sponge, &
                  "If true, sponges may be applied anywhere in the domain. \n"//&
                  "The exact location and properties of those sponges are \n"//&

src/parameterizations/vertical/MOM_diabatic_driver.F90

@@ -84,6 +84,9 @@ module MOM_diabatic_driver
   logical :: use_geothermal          !< If true, apply geothermal heating.
   logical :: use_int_tides           !< If true, use the code that advances a separate set
                                      !! of equations for the internal tide energy density.
+  logical :: ePBL_is_additive        !< If true, the diffusivity from ePBL is added to all
+                                     !! other diffusivities. Otherwise, the larger of kappa-
+                                     !! shear and ePBL diffusivities are used.
   integer :: nMode = 1               !< Number of baroclinic modes to consider
   logical :: int_tide_source_test    !< If true, apply an arbitrary generation site
                                      !! for internal tide testing (BDM)
@@ -324,6 +327,7 @@ subroutine diabatic(u, v, h, tv, fluxes, visc, ADp, CDp, dt, G, GV, CS)
 
   integer :: ig, jg      ! global indices for testing testing itide point source (BDM)
   logical :: avg_enabled ! for testing internal tides (BDM)
+  real :: Kd_add_here    ! An added diffusivity in m2/s
 
   is   = G%isc  ; ie  = G%iec  ; js  = G%jsc  ; je  = G%jec ; nz = G%ke
   Isq  = G%IscB ; Ieq = G%IecB ; Jsq = G%JscB ; Jeq = G%JecB
@@ -741,12 +745,18 @@ subroutine diabatic(u, v, h, tv, fluxes, visc, ADp, CDp, dt, G, GV, CS)
       ! Augment the diffusivities due to those diagnosed in energetic_PBL.
       do K=2,nz ; do j=js,je ; do i=is,ie
 
-        Ent_int = Kd_ePBL(i,j,K) * (GV%m_to_H**2 * dt) / &
-                  (0.5*(h(i,j,k-1) + h(i,j,k)) + h_neglect)
+        if (CS%ePBL_is_additive) then
+          Kd_add_here = Kd_ePBL(i,j,K)
+          visc%Kv_turb(i,j,K) = visc%Kv_turb(i,j,K) + Kd_ePBL(i,j,K)
+        else
+          Kd_add_here = max(Kd_ePBL(i,j,K) - visc%Kd_turb(i,j,K), 0.0)
+          visc%Kv_turb(i,j,K) = max(visc%Kv_turb(i,j,K), Kd_ePBL(i,j,K))
+        endif
+        Ent_int = Kd_add_here * (GV%m_to_H**2 * dt) / &
+                    (0.5*(h(i,j,k-1) + h(i,j,k)) + h_neglect)
         eb(i,j,k-1) = eb(i,j,k-1) + Ent_int
         ea(i,j,k) = ea(i,j,k) + Ent_int
-        visc%Kv_turb(i,j,K) = visc%Kv_turb(i,j,K) + Kd_ePBL(i,j,K)
-        Kd_int(i,j,K)  = Kd_int(i,j,K) + Kd_ePBL(i,j,K)
+        Kd_int(i,j,K)  = Kd_int(i,j,K) + Kd_add_here
 
         ! for diagnostics
         Kd_heat(i,j,K) = Kd_heat(i,j,K) + Kd_int(i,j,K)
@@ -1751,6 +1761,10 @@ subroutine diabatic_driver_init(Time, G, GV, param_file, useALEalgorithm, diag,
                  "If true, use an implied energetics planetary boundary \n"//&
                  "layer scheme to determine the diffusivity and viscosity \n"//&
                  "in the surface boundary layer.", default=.false.)
+  call get_param(param_file, mod, "EPBL_IS_ADDITIVE", CS%ePBL_is_additive, &
+                 "If true, the diffusivity from ePBL is added to all\n"//&
+                 "other diffusivities. Otherwise, the larger of kappa-\n"//&
+                 "shear and ePBL diffusivities are used.", default=.true.)
   call get_param(param_file, mod, "DOUBLE_DIFFUSION", differentialDiffusion, &
                  "If true, apply parameterization of double-diffusion.", &
                  default=.false. )

src/parameterizations/vertical/MOM_internal_tide_input.F90

@@ -325,7 +325,7 @@ subroutine int_tide_input_init(Time, G, GV, param_file, diag, CS, itide)
   CS%diag => diag
 
   ! Read all relevant parameters and write them to the model log.
-  call log_version(param_file, mod, version)
+  call log_version(param_file, mod, version, "")
 
   call get_param(param_file, mod, "INPUTDIR", CS%inputdir, default=".")
   CS%inputdir = slasher(CS%inputdir)

src/parameterizations/vertical/MOM_set_diffusivity.F90

@@ -2544,7 +2544,7 @@ subroutine set_diffusivity_init(Time, G, GV, param_file, diag, CS, diag_to_Z_CSp
 
 
   ! Read all relevant parameters and write them to the model log.
-  call log_version(param_file, mod, version)
+  call log_version(param_file, mod, version, "")
 
   call get_param(param_file, mod, "INPUTDIR", CS%inputdir, default=".")
   CS%inputdir = slasher(CS%inputdir)

src/user/user_change_diffusivity.F90

@@ -226,7 +226,7 @@ subroutine user_change_diff_init(Time, G, param_file, diag, CS)
   CS%diag => diag
 
   ! Read all relevant parameters and write them to the model log.
-  call log_version(param_file, mod, version)
+  call log_version(param_file, mod, version, "")
   call get_param(param_file, mod, "USER_KD_ADD", CS%Kd_add, &
                  "A user-specified additional diffusivity over a range of \n"//&
                  "latitude and density.", units="m2 s-1", default=0.0)