add simplified sponge layer

including namelist option to specify sponge layer bottom/start

Author: mt5555

components/eam/bld/namelist_files/namelist_defaults_eam.xml

@@ -1245,6 +1245,7 @@
 
 <!-- viscosity  -->
 
+<tom_sponge_start> 0d0 </tom_sponge_start>
 <nu_top> 2.5e5 </nu_top>
 <nu_top dyn_target="theta-l" hgrid="ne120np4"> 1e5 </nu_top>
 <nu_top dyn_target="theta-l" hgrid="ne240np4"> 4.3e4 </nu_top>

components/eam/bld/namelist_files/namelist_definition.xml

@@ -5784,6 +5784,12 @@ Second-order viscosity applied only near the model top [m^2/s].
 Default: Set by build-namelist.
 </entry>
 
+<entry id="tom_sponge_start" type="real" category="se"
+       group="ctl_nl" valid_values="" >
+Bottom of sponge layer in hPa.
+Default: 0 (use default value based on reference pressure at model top).
+</entry>
+
 <entry id="fine_ne" type="integer" category="se"
        group="ctl_nl" valid_values="" >
 baseline ne for scalar hypervis tuning

components/homme/src/share/control_mod.F90

@@ -161,6 +161,7 @@ module control_mod
   real (kind=real_kind), public :: nu_q    = -1               ! default = nu   tracer viscosity
   real (kind=real_kind), public :: nu_p    = -1               ! default = nu   ps equ. viscosity
   real (kind=real_kind), public :: nu_top  = 0.0D5            ! top-of-the-model viscosity
+  real (kind=real_kind), public :: tom_sponge_start=0         ! start of sponge layer, in hPa
 
   integer, public :: hypervis_subcycle=1                      ! number of subcycles for hyper viscsosity timestep
   integer, public :: hypervis_subcycle_tom=0                  ! number of subcycles for TOM diffusion

components/homme/src/share/namelist_mod.F90

@@ -64,6 +64,7 @@ module namelist_mod
     nu_div,        &
     nu_p,          &
     nu_top,        &
+    tom_sponge_start, &
     dcmip16_mu,     &
     dcmip16_mu_s,   &
     dcmip16_mu_q,   &
@@ -284,6 +285,7 @@ subroutine readnl(par)
       nu_div,        &
       nu_p,          &
       nu_top,        &
+      tom_sponge_start, &
       dcmip16_mu,     &
       dcmip16_mu_s,   &
       dcmip16_mu_q,   &
@@ -790,6 +792,7 @@ subroutine readnl(par)
     call MPI_bcast(nu_div,          1, MPIreal_t   , par%root,par%comm,ierr)
     call MPI_bcast(nu_p,            1, MPIreal_t   , par%root,par%comm,ierr)
     call MPI_bcast(nu_top,          1, MPIreal_t   , par%root,par%comm,ierr)
+    call MPI_bcast(tom_sponge_start,1, MPIreal_t   , par%root,par%comm,ierr)
 
     call MPI_bcast(dcmip16_mu,      1, MPIreal_t   , par%root,par%comm,ierr)
     call MPI_bcast(dcmip16_mu_s,    1, MPIreal_t   , par%root,par%comm,ierr)
@@ -1211,6 +1214,7 @@ subroutine readnl(par)
        write(iulog,'(a,2e9.2)')"viscosity:  nu_q      = ",nu_q
        write(iulog,'(a,2e9.2)')"viscosity:  nu_p      = ",nu_p
        write(iulog,'(a,2e9.2)')"viscosity:  nu_top      = ",nu_top
+       write(iulog,'(a,2e9.2)')"viscosity:  tom_sponge_start  = ",tom_sponge_start
 
        if(dcmip16_mu/=0)  write(iulog,'(a,2e9.2)')"1st order viscosity:  dcmip16_mu   = ",dcmip16_mu
        if(dcmip16_mu_s/=0)write(iulog,'(a,2e9.2)')"1st order viscosity:  dcmip16_mu_s = ",dcmip16_mu_s

components/homme/src/theta-l/share/model_init_mod.F90

@@ -23,7 +23,7 @@ module model_init_mod
   use viscosity_mod,      only: make_c0_vector
   use kinds,              only: real_kind,iulog
   use control_mod,        only: qsplit,theta_hydrostatic_mode, hv_ref_profiles, &
-       hv_theta_correction
+       hv_theta_correction, tom_sponge_start
   use time_mod,           only: timelevel_qdp, timelevel_t
   use physical_constants, only: g, TREF, Rgas, kappa
   use imex_mod,           only: test_imex_jacobian
@@ -121,35 +121,28 @@ subroutine model_init2(elem,hybrid,deriv,hvcoord,tl,nets,nete )
     nlev_tom=0
     if (hybrid%masterthread) write(iulog,*) "sponge layer nu_top viscosity scaling factor"
     do k=1,nlev
-       !press = (hvcoord%hyam(k)+hvcoord%hybm(k))*hvcoord%ps0
-       !ptop  = hvcoord%hyai(1)*hvcoord%ps0
-       ! sponge layer starts at p=4*ptop 
-       ! 
-       ! some test cases have ptop=200mb
-       if (hvcoord%etai(1)==0) then
-          ! pure sigma coordinates could have etai(1)=0
-          ptop_over_press = hvcoord%etam(1) / hvcoord%etam(k)  
+       if (tom_sponge_start==0) then
+          ! some test cases have ptop=200mb
+          if (hvcoord%etai(1)==0) then
+             ! pure sigma coordinates could have etai(1)=0
+             ptop_over_press = hvcoord%etam(1) / hvcoord%etam(k)  
+          else
+             ptop_over_press = hvcoord%etai(1) / hvcoord%etam(k)  
+          endif
+          ! active for p<10*ptop (following cd_core.F90 in CAM-FV)
+          ! CAM 26L and 30L:  top 3 levels 
+          ! E3SM 72L:  top 6 levels
+          !original cam formula
+          !nu_scale_top(k) = 8*(1+tanh(log(ptop_over_press))) ! active for p<4*ptop
+          nu_scale_top(k) = 16*ptop_over_press**2 / (ptop_over_press**2 + 1)
        else
-          ptop_over_press = hvcoord%etai(1) / hvcoord%etam(k)  
+          ptop_over_press = (tom_sponge_start/1d3) / hvcoord%etam(k)  
+          nu_scale_top(k)=0.15d0 * ptop_over_press**2
        endif
 
-       ! active for p<10*ptop (following cd_core.F90 in CAM-FV)
-       ! CAM 26L and 30L:  top 3 levels 
-       ! E3SM 72L:  top 6 levels
-       !original cam formula
-       !nu_scale_top(k) = 8*(1+tanh(log(ptop_over_press))) ! active for p<4*ptop
-       nu_scale_top(k) = 16*ptop_over_press**2 / (ptop_over_press**2 + 1)
-
        if (nu_scale_top(k)<0.15d0) nu_scale_top(k)=0
+       if (nu_scale_top(k)>8d0) nu_scale_top(k)=8d0
 
-       !nu_scale_top(k) = 8*(1+.911*tanh(log(ptop_over_press))) ! active for p<6.5*ptop
-       !if (nu_scale_top(k)<1d0) nu_scale_top(k)=0
-
-       ! original CAM3/preqx formula
-       !if (k==1) nu_scale_top(k)=4
-       !if (k==2) nu_scale_top(k)=2
-       !if (k==3) nu_scale_top(k)=1
-       !if (k>3) nu_scale_top(k)=0
 
        if (nu_scale_top(k)>0) nlev_tom=k