anchor.diff:  
               ! Still not 100% convinced this works well in cases other than flat_z_height = 0 (w sleve). So for now best to keep at 0 when using sleve?
               max_level = find_flat_model_level(options, nz, dz)
  -            
  +
               if(max_level /= nz) then
                   if (this_image()==1) then
                       print*, "    flat z height ", options%parameters%flat_z_height
                       print*, "    flat z height set to 0 to comply with SLEVE coordinate calculation "
                       print*, "    flat z height now", nz
                   end if
  -                max_level = nz  
  +                max_level = nz
               end if
   
               smooth_height = sum(dz(1:max_level)) !sum(global_terrain) / size(global_terrain) + sum(dz(1:max_level))
  
positive.code:             s1 = smooth_height / options%parameters%decay_rate_L_topo
              s2 = smooth_height / options%parameters%decay_rate_S_topo
              n  =  options%parameters%sleve_n  ! this will have an effect on the z_level ratio throughout the vertical column, and thus on the terrain induced acceleration with wind=2 . Conceptually very nice, but for wind is 2 not ideal. Unless we let that acceleration depend on the difference between hi-res and lo-res terrain.
              
  
              ! Scale dz with smooth_height/sum(dz(1:max_level)) before calculating sleve levels.
              dz_scl(:)   =   dz(1:nz) ! *  smooth_height / sum(dz(1:max_level))  ! this leads to a jump in dz thickness at max_level+1. Not sure if this is a problem.
              
  
              ! - - -   calculate invertibility parameter gamma (Schär et al 2002 eqn 20):  - - - - - -
              gamma  =  1  -  MAXVAL(h1)/s1 * COSH(smooth_height/s1)/SINH(smooth_height/s1) - MAXVAL(h2)/s2 * COSH(smooth_height/s2)/SINH(smooth_height/s2)
  
              ! COSMO1 operational setting (but model top is at ~22000 masl):
              !    Decay Rate for Large-Scale Topography: svc1 = 10000.0000  
              !    Decay Rate for Small-Scale Topography: svc2 =  3300.0000
              if ((this_image()==1)) then
                  print*, "    Using a SLEVE coordinate with a Decay height for Large-Scale Topography: (s1) of ", s1, " m."
  
negative.code: 
  
      !> -------------------------------
      !! Setup the SLEVE vertical grid structure.
      !!   This basically entails 2 transformations: First a linear one so that sum(dz) ranges from 0 to smooth_height H.
      !!   (boundary cnd (3) in Schär et al 2002)  Next, the nonlinear SLEVE transformation
      !!    eqn (2) from Leuenberger et al 2009 z_sleve = Z + terrain * sinh((H/s)**n - (Z/s)**n) / SINH((H/s)**n) (for both smallscale and largescale terrain)
  