Add changes for variable epssm and reference sounding option

[jklemp]

1. Code added to ../src/core_atmosphere/mpas_atm_core.F to compute
   vertical profiles of epssm that vary with height at the start of the
   time integration. These profiles are determined by the parameters
   config_min_coeff, config_max_coeff, config_transition_lower_bound,
   and config_transition_upper_bound, and are specified in namelist.atmosphere.

2. Code added to allow thermodynamic variables to be represented either
   as perurbations from a specified reference state or as the full
   variables without a reference state. This option is controlled by the
   logical config_reference_sounding as specified in namelist.init_atmosphere.
   When configured using the full variables, the horizontal and vertical
   pressure gradients computed on the large time steps are written in terms
   of log(p).

[skamaroc]

In src/core_atmosphere/Registry.xml

config_min_coeff -> config_epssm_minimum
config_max_coeff -> config_epssm_maximum
config_transition_lower_bound -> config_epssm_transition_bottom_z
units=“-“ -> units=“m”
config_transition_upper_bound -> config_epssm_transition_top_z
units=“-“ -> units=“m”

(Also - code changes to handle the namelist variable name changes)

Change descriptions in Registry
<var name="etp" type="real" dimensions="nVertLevels" units="unitless"
description="etp = (1+epssm(z))/2 at theta point”/>
<var name="etm” type="real" dimensions="nVertLevels" units="unitless"
description="etm = (1-epssm(z))/2 at theta point”/>
<var name="ewp" type="real" dimensions="nVertLevelsP1” units="unitless"
description="ewp = (1+epssm(z))/2 at w point”/>
<var name="ewm” type="real" dimensions="nVertLevelsP1” units="unitless"

•        description="ewm = (1-epssm(z))/2 at w point”/>
  

In src/core_atmosphere/mpas_atm_core.F

Let’s delete the following commented lines

! min_coeff = 0.1_RKIND
! max_coeff = 0.5_RKIND
! transition_lower_bound = 8000._RKIND
! transition_upper_bound = 13000._RKIND

! call mpas_log_write('test $i $r', intARgs=(/k/),realArgs=(/rdzw(k)/))

In src/core_atmosphere/dynamics/mpas_atm_time_integration.F

Let’s delete my (WCS_epssm) old code in subroutine atm_srk3:

!----------------------
!SK_epssm: setting the 'mesh' pool pointer
call mpas_pool_get_subpool(domain % blocklist % structs, 'mesh', mesh)
!----------------------------

!------------------------
!SK_epssm: added this to retrieve the pointer for rdzw
call mpas_pool_get_array(mesh, 'rdzw', rdzw)
!------------------

  ! WCS_epssm                                                                                                                               
  ! allocate storage for epssm column arrays and set values                                                                                 

! allocate(ewp(nVertLevels+1))
! allocate(ewm(nVertLevels+1))
! allocate(etp(nVertLevels))
! allocate(etm(nVertLevels))
! call set_epssm(ewp,ewm,etp,etm,nVertLevels,rdzw)

  ! WCS_epssm                                                                                                                               

! deallocate(ewp)
! deallocate(ewm)
! deallocate(etp)
! deallocate(etm)

In subroutine atm_compute_vert_imp_coefs_work

Let’s get rid of the test print

! temporary print to check etm and ewm

  integer :: iprint=0      
  if(iprint.eq.0) then
     do k=1,nVertLevels
        call mpas_log_write(' k, etp(k), ewp(k) = $i $r $r', intArgs=(/k/), realArgs=(/etp(k),ewp(k)/)) 
     end do
     iprint = 1
  end if 

Let’s get rid of old WCS code that commented out

! WCS_epssm
! cofwr(k,iCell) =.5dtsepsgravity*(fzm(k)*zz(k,iCell)+fzp(k)*zz(k-1,iCell))

        ! WCS_epssm                                                                                                                         
        ! cofwz(k,iCell) = dtseps*c2*(fzm(k)*zz(k,iCell)+fzp(k)*zz(k-1,iCell))  &                                                           
        !      *rdzu(k)*cqw(k,iCell)*(fzm(k)*p (k,iCell)+fzp(k)*p (k-1,iCell))                                                              
        ! coftz(k,iCell) = dtseps*   (fzm(k)*t (k,iCell)+fzp(k)*t (k-1,iCell)) 


        ! WCS_epssm                                                                                                                         
        ! cofwt(k,iCell) = .5*dtseps*rcv*zz(k,iCell)*gravity*rb(k,iCell)/(1.+qtotal)  &                                                     
        !                     *p(k,iCell)/((rtb(k,iCell)+rt(k,iCell))*pb(k,iCell))   


        ! WiCS_epssm                                                                                                                        
        ! b_tri(k) = 1.                                                  &                                                                  
        !              +cofwz(k  ,iCell)*(coftz(k  ,iCell)*rdzw(k  )*zz(k  ,iCell)   &                                                      
        !                           +coftz(k  ,iCell)*rdzw(k-1)*zz(k-1,iCell))   &                                                          
        !              -coftz(k  ,iCell)*(cofwt(k  ,iCell)*rdzw(k  )             &                                                          
        !                            -cofwt(k-1,iCell)*rdzw(k-1))            &                                                              
        !              +cofwr(k,  iCell)*(cofrz(k    )-cofrz(k-1))              


        ! WCS_epssm                                                                                                                         
        ! alpha_tri(k,iCell) = 1./(b_tri(k)-a_tri(k,iCell)*gamma_tri(k-1,iCell))                                                            
        ! gamma_tri(k,iCell) = c_tri(k)*alpha_tri(k,iCell)  

In subroutine atm_advance_acoustic_step_work

Let’s get rid of the print statement

! temporary print to check etm and ewm

  integer :: iprint=0      
  if(iprint.eq.0) then
     do k=1,nVertLevels
        call mpas_log_write(' k, etm(k), ewm(k) = $i $r $r', intArgs=(/k/), realArgs=(/etm(k),ewm(k)/)) 
     end do
     iprint = 1
  end if 

Let’s get rid of the old commented out WCS code

        ! WCS_epssm                                                                                                                         
        ! rs(k) = rho_pp(k,iCell) + dts*tend_rho(k,iCell) + rs(k)                  &                                                        
        !                 - cofrz(k)*resm*(rw_p(k+1,iCell)-rw_p(k,iCell))          

        ! WCS_epssm                                                                                                                         
        ! ts(k) = rtheta_pp(k,iCell) + dts*tend_rt(k,iCell) + ts(k)                &                                                        
        !                    - resm*rdzw(k)*( coftz(k+1,iCell)*rw_p(k+1,iCell)     &                                                        
        !                                    -coftz(k,iCell)*rw_p(k,iCell))        

        ! WCS_epssm                                                                                                                         
        ! wwavg(k,iCell) = wwavg(k,iCell) + 0.5*(1.0-epssm)*rw_p(k,iCell)  

       ! WCS_epssm                                                                                                                          
        ! rw_p(k,iCell) = rw_p(k,iCell) +  dts*tend_rw(k,iCell)                       &                                                     
        !            - cofwz(k,iCell)*((zz(k  ,iCell)*ts(k)                           &                                                     
        !                          -zz(k-1,iCell)*ts(k-1))                            &                                                     
        !                    +resm*(zz(k  ,iCell)*rtheta_pp(k  ,iCell)                &                                                     
        !                          -zz(k-1,iCell)*rtheta_pp(k-1,iCell)))              &                                                     
        !            - cofwr(k,iCell)*((rs(k)+rs(k-1))                                &                                                     
        !                    +resm*(rho_pp(k,iCell)+rho_pp(k-1,iCell)))               &                                                     
        !            + cofwt(k  ,iCell)*(ts(k  )+resm*rtheta_pp(k  ,iCell))           &                                                     
        !            + cofwt(k-1,iCell)*(ts(k-1)+resm*rtheta_pp(k-1,iCell))         

       ! WCS_epssm                                                                                                                          
        ! rw_p(k,iCell) = (rw_p(k,iCell)-a_tri(k,iCell)*rw_p(k-1,iCell))*alpha_tri(k,iCell)  

       ! WCS_epssm                                                                                                                          
        ! wwAvg(k,iCell) = wwAvg(k,iCell) + 0.5*(1.0+epssm)*rw_p(k,iCell) 

       ! WCS_epssm                                                                                                                          
        ! rho_pp(k,iCell) = rs(k) - cofrz(k) *(rw_p(k+1,iCell)-rw_p(k  ,iCell))                                                             
        ! rtheta_pp(k,iCell) = ts(k) - rdzw(k)*(coftz(k+1,iCell)*rw_p(k+1,iCell)  &                                                         
        !                    -coftz(k  ,iCell)*rw_p(k  ,iCell))   

           ! WCS_epssm                                                                                                                      
           ! wwAvg(k,iCell) = wwAvg(k,iCell) + 0.5*(1.0+epssm)*rw_p(k,iCell)      

Let’s get rid of subroutine set_epssm(ewp,ewm,etp,etm,n,rdzw)
And subroutine calculate_computational_height(nVertLevels,rdzw,height_u_levels, height_w_levels)