Merge branch 'mt5555/tom_sponge_start' (PR #4723)

Add a simplified sponge layer with namelist variable (tom_sponge_start) to set the starting layer

This is needed to run at high resolution with the E3SM 72L configuration (72L sponge layer is too thin, it starts at the 1mb layer). 72L configurations will crash at NE256 and higher resolution without making the sponge layer thicker.

It will also be useful for the SCREAM 128L configuration, where the sponge layer is most likely too thick (it starts at the 20mb layer)

default: tom_sponge_start=0 (disabled).

documentation:
https://acme-climate.atlassian.net/wiki/spaces/DOC/pages/2967798203/EAM+Top+of+Model+Sponge+Layer

[BFB]

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,   &
@@ -285,6 +286,7 @@ subroutine readnl(par)
       nu_div,        &
       nu_p,          &
       nu_top,        &
+      tom_sponge_start, &
       dcmip16_mu,     &
       dcmip16_mu_s,   &
       dcmip16_mu_q,   &
@@ -793,6 +795,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)
@@ -1208,6 +1211,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
@@ -119,35 +119,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