anchor.diff:  
               ! Offset analogous to: z_u(:,i,:) = z_u(:,i,:) + dz(i) / 2 * zr_u(:,i,:)
               z_u(:,i,:)  = dz_scl(i)/2  &
  -                          ! + terrain_u  *  SINH( (H/s)**n - ( dz_scl(i)/2 /s)**n ) / SINH((H/s)**n)
  -                          + h1_u  *  SINH( (smooth_height/s1)**n - (dz_scl(i)/2/s1)**n ) / SINH((smooth_height/s1)**n)  &! large-scale terrain
  -                          + h2_u  *  SINH( (smooth_height/s2)**n - (dz_scl(i)/2/s2)**n ) / SINH((smooth_height/s2)**n)   ! small terrain features
  +                        + h1_u  *  SINH( (smooth_height/s1)**n - (dz_scl(i)/2/s1)**n ) / SINH((smooth_height/s1)**n)  &! large-scale terrain
  +                        + h2_u  *  SINH( (smooth_height/s2)**n - (dz_scl(i)/2/s2)**n ) / SINH((smooth_height/s2)**n)   ! small terrain features
               z_v(:,i,:)  = dz_scl(i)/2   &
  -                          ! + terrain_v  *  SINH( (H/s)**n - ( dz_scl(i)/2 /s)**n ) / SINH((H/s)**n)
  -                          + h1_v  *  SINH( (smooth_height/s1)**n - (dz_scl(i)/2/s1)**n ) / SINH((smooth_height/s1)**n)  &! large-scale terrain
  -                          + h2_v  *  SINH( (smooth_height/s2)**n - (dz_scl(i)/2/s2)**n ) / SINH((smooth_height/s2)**n)   ! small terrain features
  +                        + h1_v  *  SINH( (smooth_height/s1)**n - (dz_scl(i)/2/s1)**n ) / SINH((smooth_height/s1)**n)  &! large-scale terrain
  +                        + h2_v  *  SINH( (smooth_height/s2)**n - (dz_scl(i)/2/s2)**n ) / SINH((smooth_height/s2)**n)   ! small terrain features
   
               zr_u(:,i,:)  =  (z_u(:,i,:) - terrain_u) / ( dz_scl(i)/2 )
               zr_v(:,i,:)  =  (z_v(:,i,:) - terrain_v) / (dz_scl(i)/2 )
  
positive.code:             z_u(:,i,:)          = z_u(:,i,:) + dz(i) / 2 * zr_u(:,i,:)
              z_v(:,i,:)          = z_v(:,i,:) + dz(i) / 2 * zr_v(:,i,:)
  
              do i = this%grid%kms+1, this%grid%kme
                  if (i<=max_level) then
                      jacobian(:,i,:) = jacobian(:,i-1,:)
  
negative.code: 
              i = this%grid%kme + 1
              global_z_interface(:,i,:) = global_z_interface(:,i-1,:) + global_dz_interface(:,i-1,:)
          end associate
          
      end subroutine setup_simple_z
  
  
  
      !> -------------------------------
      !! Initialize various domain variables, mostly z, dz, etc.
      !!
      !! -------------------------------
      subroutine initialize_core_variables(this, options)
          implicit none
          class(domain_t), intent(inout)  :: this
          type(options_t), intent(in)     :: options
  
          real, allocatable :: temp(:,:,:)
  
          call read_core_variables(this, options)
  
          call allocate_z_arrays(this)
  
          ! Setup the vertical grid structure, either as a SLEVE coordinate, or a more 'simple' vertical structure:
          if (options%parameters%sleve) then 
              
              call split_topography(this, options)  ! here h1 and h2 are calculated
              call setup_sleve(this, options)
         
          else  
         
              call setup_simple_z(this, options)
          
          endif 
  
          associate(ims => this%ims,      ime => this%ime,                        &
                    jms => this%jms,      jme => this%jme,                        &
                    kms => this%kms,      kme => this%kme,                        &
                    z                     => this%z%data_3d,                      &
                    global_jacobian       => this%global_jacobian,                &
                    jacobian              => this%jacobian,                       &
                    jacobian_u            => this%jacobian_u,                     &
                    jacobian_v            => this%jacobian_v,                     &
                    jacobian_w            => this%jacobian_w,                     &
                    zr_u                  => this%zr_u,                           &
                    zr_v                  => this%zr_v)
  
  
              if (allocated(temp)) deallocate(temp)
              allocate(temp(this%ids:this%ide+1, this%kds:this%kde, this%jds:this%jde+1))
              temp(this%ids,:,this%jds:this%jde) = global_jacobian(this%ids,:,this%jds:this%jde)
              temp(this%ide+1,:,this%jds:this%jde) = global_jacobian(this%ide,:,this%jds:this%jde)
              temp(this%ids+1:this%ide,:,this%jds:this%jde) = (global_jacobian(this%ids+1:this%ide,:,this%jds:this%jde) + &
                                                                  global_jacobian(this%ids:this%ide-1,:,this%jds:this%jde))/2
              jacobian_u = temp(ims:ime+1,:,jms:jme)
  
              temp(this%ids:this%ide,:,this%jds) = global_jacobian(this%ids:this%ide,:,this%jds)
              temp(this%ids:this%ide,:,this%jde+1) = global_jacobian(this%ids:this%ide,:,this%jde)
              temp(this%ids:this%ide,:,this%jds+1:this%jde) = (global_jacobian(this%ids:this%ide,:,this%jds+1:this%jde) + &
                                                  global_jacobian(this%ids:this%ide,:,this%jds:this%jde-1))/2
              jacobian_v = temp(ims:ime,:,jms:jme+1)
  
              temp(this%ids:this%ide,this%kme,this%jds) = global_jacobian(this%ids:this%ide,this%kme,this%jds)
              temp(this%ids:this%ide,this%kms:this%kme-1,this%jds:this%jde) = (global_jacobian(this%ids:this%ide,this%kms:this%kme-1,this%jds:this%jde) + &
                                                                              global_jacobian(this%ids:this%ide,this%kms+1:this%kme,this%jds:this%jde))/2
              jacobian_w = temp(ims:ime,:,jms:jme)
  
              call setup_dzdxy(this, options)
  
                  ! technically these should probably be defined to the k+1 model top as well bu not used at present.
                  ! z_interface(:,i,:) = z_interface(:,i-1,:) + dz_interface(:,i-1,:)
          end associate
  
          ! z_u and zr_u are on the v/u_grid2d_ext; move to vu_grid2d
  