anchor.diff:              call setup_sleve(this, options)
   
           else
  -
  +            ! This will set up either a Gal-Chen terrainfollowing coordinate, or no terrain following. 
               call setup_simple_z(this, options)
   
           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
  +        ! 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,                        &
  
positive.code:         h2_v =  h_v  - h1_v
  
          ! In case one wants to see how the terrain is split by smoothing, activate the block below and run in debug:
          if ((this_image()==1).and.(options%parameters%debug)) then
            call io_write("terrain_smooth_h1.nc", "h1", h1(:,:) )
            call io_write("terrain_smooth_h2.nc", "h2", h2(:,:) )
          !   call io_write("h1_u.nc", "h1_u", h1_u(:,:) )
          !   call io_write("h2_u.nc", "h2_u", h2_u(:,:) )
          endif
  
          if (this_image()==1) then
             print*, "       Max of full topography", MAXVAL(global_terrain )
  
negative.code:                 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_min
                  if(gamma_min <= 0) print*, " CAUTION: coordinate transformation is not invertible (gamma <= 0 ) !!! reduce decay rate(s), and/or increase flat_z_height!"
                  ! if(options%parameters%debug)  write(*,*) "   (for (debugging) reference: 'gamma(n=1)'= ", gamma,")"
                  print*, ""
              endif
  
  