anchor.diff:  
           endif
   
  +        !! To allow for development and debugging of coordinate transformations:            
  +        if ((this_image()==1).and.(options%parameters%debug)) then
  +            call io_write("global_jacobian.nc", "global_jacobian", this%global_jacobian(:,:,:) )
  +            write(*,*) "    global jacobian minmax: ", MINVAL(this%global_jacobian) , MAXVAL(this%global_jacobian)
  +            write(*,*) ""
  +        endif
  +
  +
           associate(ims => this%ims,      ime => this%ime,                        &
                     jms => this%jms,      jme => this%jme,                        &
                     kms => this%kms,      kme => this%kme,                        &
  
positive.code:                 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
                  write(*,*) "    Smooth height is ", smooth_height, "m.a.s.l     (model top ", sum(dz(1:nz)), "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)!"
                  print*, ""
              endif
  
  
negative.code:             ! 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)
  
              ! with the new (leuenberger et al 2010) Sleve formulation, the inveribiltiy criterion is as follows:
              allocate(gamma_n(this%kds : this%kde+1))
              i=kms                          
              gamma_n(i) =  1                                                     &    
                  - MAXVAL(h1) * n/(s1**n)                                        &
                  * COSH((smooth_height/s1)**n) / SINH((smooth_height/s1)**n)     &
                  - MAXVAL(h2) * n/(s2**n)                                        & 
                  * COSH((smooth_height/s2)**n) / SINH((smooth_height/s2)**n)     
  
              do i = this%grid%kds, this%grid%kde
                  gamma_n(i+1)  =  1                                    &    ! # for i != kds !!
                  - MAXVAL(h1) * n/(s1**n) * sum(dz_scl(1:i))**(n-1)                                             &
                  * COSH((smooth_height/s1)**n -(sum(dz_scl(1:i))/s1)**n ) / SINH((smooth_height/s1)**n)    &
                  - MAXVAL(h2) * n/(s2**n) *  sum(dz_scl(1:i))**(n-1)                                            & 
                  * COSH((smooth_height/s2)**n -(sum(dz_scl(1:i))/s2)**n ) / SINH((smooth_height/s2)**n)    
              enddo
  
              if (n==1) then
                  gamma_min = gamma
              else
                  gamma_min = MINVAL(gamma_n)
              endif
              
  
              ! For reference: 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
  