anchor.diff:              z_u(:,i,:)          = z_u(:,i,:) + dz(i) / 2 * zr_u(:,i,:)
               z_v(:,i,:)          = z_v(:,i,:) + dz(i) / 2 * zr_v(:,i,:)
   
  +            ! Now the higher (k!=1) levels can be calculated:
               do i = this%grid%kms+1, this%grid%kme
                   if (i<=max_level) then
                       jacobian(:,i,:) = jacobian(:,i-1,:)
  
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.
              ! h = 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.
              ! dz_scl(:)   =   dz(1:nz)  *  H / sum(dz(1:nz))  ! gives the same for flatz=0, but smoother otherwise? BAD idea
              ! dz_scl   =   dz(:)  *  smooth_height / sum(dz(1:max_level))
              ! 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(smooth_height/s1)/SINH(smooth_height/s1) - MAXVAL(h2)/s2 * COSH(smooth_height/s2)/SINH(smooth_height/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
                ! print*, "using a sleve_decay_factor (H/s) of ", options%parameters%sleve_decay_factor
                print*, "    Using a SLEVE coordinate with a Decay height for Large-Scale Topography: (s1) of ", s1, " m."
                print*, "    Using a SLEVE coordinate with a Decay height for Small-Scale Topography: (s2) of ", s2, " m."
                print*, "    Using a sleve_n of ", options%parameters%sleve_n
                ! print*, ""
                write(*,*) "    Smooth height (model top) is ", smooth_height, "m.a.s.l"
                write(*,*) "    invertibility parameter gamma is: ", gamma
                if(gamma <= 0) print*, " CAUTION: coordinate transformation is not invertible (gamma <= 0 ) !!! reduce decay rate(s)!"
                ! write(*,*) "  mean terrain ", sum(terrain) / size(terrain)
                ! write(*,*) "  sum(dz) ", sum(dz(1:max_level))
                ! write(*,*) "  sum(dz_scl) ", sum(dz_scl(1:max_level))
                ! write(*,*) "  model top ", sum(dz_scl(1:nz))
                print*, ""
  
              endif
  
              i=kms
  
              ! - - - - -   Mass grid calculations for lowest level (i=kms)  - - - - -
  
              !use temp to store global z-interface so that global-jacobian can be calculated
              allocate(temp(this%ids:this%ide, this%kds:this%kde, this%jds:this%jde))
  
              temp(:,i,:)   = global_terrain
  
negative.code:         endif
  
  
          ! create a separate variable that will be smoothed later on:
          h1 =  global_terrain(this%grid2d%ids:this%grid2d%ide, this%grid2d%jds:this%grid2d%jde)
  
          ! offset the global terrain for the h_(u/v) calculations:
          call array_offset_x_2d(global_terrain, temp_offset)
          h_u = temp_offset
          h1_u = temp_offset
          if (allocated(temp_offset)) deallocate(temp_offset)
  
          call array_offset_y_2d(global_terrain, temp_offset)
          h_v = temp_offset
          h1_v = temp_offset
  
  