anchor.diff:              ! 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
  
positive.code:                   jms => this%jms,      jme => this%jme,                        &
                    kms => this%kms,      kme => this%kme,                        &
                    terrain               => this%terrain%data_2d,                &
                    terrain_u             => this%terrain_u,                      &
                    terrain_v             => this%terrain_v,                      &
                    forcing_terrain       => this%forcing_terrain%data_2d,        &
  
negative.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 ! only for the -currently unused- delta_dzdx_sc calculation (1B)
          ! wind_top = (s1+s2)/2 ! Experiment, lets see what this does. 
  
          ! To prevent the wind_top (the height below which we hor.accelerate winds) from becoming too low, thus creating 
          !   very large acceleration, we introduce this check. 
          e = 1.2  ! <- first guess
          if (MAXVAL(global_terrain) *e < s1 ) then
              wind_top = s1 
              if (this_image()==1) print*, "  horizontally accelerating winds below:", wind_top, "m. " !,"(Factor H/s:", H/s ,")"
          else
              wind_top = MAXVAL(global_terrain) * e !**2
              if (this_image()==1 )   print*, "  adjusting wind top upward from ",s1 ," to ", wind_top  ,"m. Horizontally accelerating winds below this level."
          endif
          ! if (this_image()==1) print*, "  s_accel max: ", wind_top, "  - h max:", MAXVAL(global_terrain)
  
  
          !_________ 1. Calculate delta_dzdx for w_real calculation - CURRENTLY NOT USED- reconsider  _________
          allocate(delta_terrain(this% ims : this% ime, &
  