anchor.diff:  
           real, allocatable :: temp(:,:,:) !, terrain_u(:,:), terrain_v(:,:)
           integer :: i
  -        real :: smooth_height, H, s, n, s1, s2
  +        real :: smooth_height, H, s, n, s1, s2, gamma
           logical :: SLEVE  
           ! character :: filename, file_idS, file_idn
   
  
positive.code:         ! s  =  H / options%parameters%sleve_decay_factor  
          s1 =  H / options%parameters%decay_rate_L_topo 
          s2 =  H / options%parameters%decay_rate_S_topo  
  
          s = (s1+s2)/2 ! Experiment, lets see what this does. 
          ! if (this_image()==1) print*, "  s_accel max: ", s, "  - h max:", MAXVAL(global_terrain)
            
  
          !_________ 1. Calculate delta_dzdx for w_real calculation - CURRENTLY NOT USED- reconsider  _________
          allocate(delta_terrain(this% ims : this% ime, &
  
negative.code:             ! H        =  sum(dz_scl(1:max_level))  ! should also lead to smooth_height, but more error proof?
  
  
              ! - - -   calculate invertibility parameter gamma (Schär et al 2002 eqn 20):  - - - - - -
              gamma  =  1  -  MAXVAL(h1)/s1 * COSH(H/s1)/SINH(H/s1) - MAXVAL(h2)/s2 * COSH(H/s2)/SINH(H/s2)
  
  
              ! Decay Rate for Large-Scale Topography: svc1 = 10000.0000  COSMO1 operational setting (but model top is at ~22000 masl)
              ! Decay Rate for Small-Scale Topography: svc2 =  3300.0000
              if ((this_image()==1)) then
  