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:         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:                     zr_u(:,i,:)  = (z_u(:,i,:) - z_u(:,i-1,:)) / (dz_scl(i)/2 + dz_scl(i-1)/2 )  ! if dz_scl(i-1) = 0 (and no error)  k=1 can be included
                      zr_v(:,i,:)  = (z_v(:,i,:) - z_v(:,i-1,:)) / (dz_scl(i)/2 + dz_scl(i-1)/2 )
  
                      ! TODO does this need to go the parallelized grid (i.e ims:ime like line 1060??)
  
                  else ! above the flat_z_height
  
  