Merge branch 'atmosphere/relax_zone_divdamp_coef' into develop (PR #1082)

This merge enables namelist control of the coefficient used in the Laplacian
filter for momentum in the relaxation zone.

The Laplacian filter for momentum in the relaxation zone previously scaled the
divergent component of the Laplacian with a hard-wired coefficient. With this
merge, a new namelist option, config_relax_zone_divdamp_coef, has been
introduced to allow for namelist control over this coefficient. The default
value of config_relax_zone_divdamp_coef is 6.0, which is different from the
previously hard-wired coefficient of 4.0 in the
atm_bdy_adjust_dynamics_relaxzone_tend routine, reflecting the outcome of more
recent testing.

Author: mgduda

src/core_atmosphere/Registry.xml

@@ -246,6 +246,11 @@
                      units="-"
                      description="Number of halo layers for fields"
                      possible_values="Integer values, typically 2 or 3; DO NOT CHANGE"/>
+
+                <nml_option name="config_relax_zone_divdamp_coef" type="real" default_value="6.0" in_defaults="false"
+                     units="-"
+                     description="Coefficient for the divergent component of the Laplacian filter of momentum in the relaxation zone"
+                     possible_values="Positive real values"/>
         </nml_record>
 
         <nml_record name="damping" in_defaults="true">

src/core_atmosphere/dynamics/mpas_atm_time_integration.F

@@ -682,7 +682,7 @@ subroutine atm_srk3(domain, dt, itimestep, exchange_halo_group)
 
 !$OMP PARALLEL DO
                do thread=1,nThreads
-                  call atm_bdy_adjust_dynamics_relaxzone_tend( tend, state, diag, mesh, nVertLevels, dt,                   &
+                  call atm_bdy_adjust_dynamics_relaxzone_tend( block % configs, tend, state, diag, mesh, nVertLevels, dt,  &
                                                                ru_driving_values, rt_driving_values, rho_driving_values,   &
                                                                cellThreadStart(thread), cellThreadEnd(thread),             &
                                                                edgeThreadStart(thread), edgeThreadEnd(thread),             &
@@ -5761,7 +5761,7 @@ end subroutine atm_bdy_adjust_dynamics_speczone_tend
 
 !-------------------------------------------------------------------------
 
-   subroutine atm_bdy_adjust_dynamics_relaxzone_tend( tend, state, diag, mesh, nVertLevels, dt,                                             &
+   subroutine atm_bdy_adjust_dynamics_relaxzone_tend( config, tend, state, diag, mesh, nVertLevels, dt,                                             &
                                                       ru_driving_values, rt_driving_values, rho_driving_values,                                     &
                                                       cellStart, cellEnd, edgeStart, edgeEnd,                               &
                                                       cellSolveStart, cellSolveEnd, edgeSolveStart, edgeSolveEnd )
@@ -5774,6 +5774,7 @@ subroutine atm_bdy_adjust_dynamics_relaxzone_tend( tend, state, diag, mesh, nVer
       !
       !  WCS Fall 2016
 
+      type (mpas_pool_type), intent(in) :: config
       type (mpas_pool_type), intent(in) :: state
       type (mpas_pool_type), intent(inout) :: tend
       type (mpas_pool_type), intent(in) :: diag
@@ -5795,6 +5796,7 @@ subroutine atm_bdy_adjust_dynamics_relaxzone_tend( tend, state, diag, mesh, nVer
 
 
       real (kind=RKIND) :: edge_sign, laplacian_filter_coef, rayleigh_damping_coef, r_dc, r_dv, invArea
+      real (kind=RKIND), pointer :: divdamp_coef
       real (kind=RKIND), dimension(nVertLevels) :: divergence1, divergence2, vorticity1, vorticity2      
       integer :: iCell, iEdge, i, k, cell1, cell2, iEdge_vort, iEdge_div
       integer :: vertex1, vertex2, iVertex
@@ -5829,6 +5831,8 @@ subroutine atm_bdy_adjust_dynamics_relaxzone_tend( tend, state, diag, mesh, nVer
       call mpas_pool_get_array(mesh, 'edgesOnVertex', edgesOnVertex)
       call mpas_pool_get_array(mesh, 'nEdgesOnCell',nEdgesOnCell)
       call mpas_pool_get_array(mesh, 'verticesOnEdge', verticesOnEdge)
+
+      call mpas_pool_get_config(config, 'config_relax_zone_divdamp_coef', divdamp_coef)
 
       !  First, Rayleigh damping terms for ru, rtheta_m and rho_zz
 
@@ -5940,8 +5944,9 @@ subroutine atm_bdy_adjust_dynamics_relaxzone_tend( tend, state, diag, mesh, nVer
             ! Compute diffusion, computed as \nabla divergence - k \times \nabla vorticity
             !
              do k=1,nVertLevels
-                tend_ru(k,iEdge) = tend_ru(k,iEdge) + laplacian_filter_coef * (  4.0*( divergence2(k) - divergence1(k) ) * r_dc  &
-                                                                                    -( vorticity2(k)  - vorticity1(k)  ) * r_dv )
+                tend_ru(k,iEdge) = tend_ru(k,iEdge) + laplacian_filter_coef &
+                                                      * (divdamp_coef * (divergence2(k) - divergence1(k)) * r_dc  &
+                                                                       -(vorticity2(k)  - vorticity1(k))  * r_dv)
              end do
 
           end if  !  end test for relaxation-zone edge