anchor.diff:  
                       if (i==this%grid%kme) then  ! if we are at the model top i+1 is not defined
   
  -                      dz_interface(:,i,:)  =  smooth_height - z_interface(:,i,:)
  -                      global_dz_interface(:,i,:)  =  smooth_height - temp(:,i,:)
  +                    dz_interface(:,i,:)  =  smooth_height - z_interface(:,i,:)
  +                    global_dz_interface(:,i,:)  =  smooth_height - temp(:,i,:)
                       else
   
  -                      temp(:,i+1,:)  = sum(dz_scl(1:i))   &
  +                    temp(:,i+1,:)  = sum(dz_scl(1:i))   &
                                       + h1  *  SINH( (smooth_height/s1)**n - (sum(dz_scl(1:i))/s1)**n ) / SINH((smooth_height/s1)**n)  &! large-scale terrain
                                       + h2  *  SINH( (smooth_height/s2)**n - (sum(dz_scl(1:i))/s2)**n ) / SINH((smooth_height/s2)**n)   ! small terrain features
  -                                     ! + terrain  *  SINH( (H/s)**n - (dz_scl(i)/s)**n ) / SINH((H/s)**n)
  -                      z_interface(:,i+1,:) = temp(ims:ime,i+1,jms:jme)
  -
  -                      global_dz_interface(:,i,:)  =  temp(:,i+1,:) - temp(:,i,:)
  -                      dz_interface(:,i,:)  =  z_interface(:,i+1,:) - z_interface(:,i,:)
  +                                    
  +                    z_interface(:,i+1,:) = temp(ims:ime,i+1,jms:jme)
  +
  +                    global_dz_interface(:,i,:)  =  temp(:,i+1,:) - temp(:,i,:)
  +                    dz_interface(:,i,:)  =  z_interface(:,i+1,:) - z_interface(:,i,:)
   
                       endif
   
                       if ( ANY(dz_interface(:,i,:)<0) ) then   ! Eror catching. Probably good to engage.
  -                      if (this_image()==1) then
  +                    if (this_image()==1) then
                           write(*,*) "Error: dz_interface below zero (for level  ",i,")"
                           print*, "min max dz_interface: ",MINVAL(dz_interface(:,i,:)),MAXVAL(dz_interface(:,i,:))
                           error stop
                           print*, dz_interface(:,i,:)
                           print*,""
  -                      endif
  +                    endif
                       else if ( ANY(dz_interface(:,i,:)<=0.01) ) then
  -                      if (this_image()==1)  write(*,*) "WARNING: dz_interface very low (at level ",i,")"
  +                    if (this_image()==1)  write(*,*) "WARNING: dz_interface very low (at level ",i,")"
                       endif
   
                       ! - - - - -   u/v grid calculations - - - - -
                       z_u(:,i,:)   = (sum(dz_scl(1:(i-1))) + dz_scl(i)/2)   &
  -                                   ! + terrain_u  *  SINH( (H/s)**n - ( (sum(dz_scl(1:(i-1)))+dz_scl(i)/2) /s)**n ) / SINH((H/s)**n)
  -                                   + h1_u  *  SINH( (smooth_height/s1)**n -  ( (sum(dz_scl(1:(i-1)))+dz_scl(i)/2) /s1)**n ) / SINH((smooth_height/s1)**n)  &! large-scale terrain
  -                                   + h2_u  *  SINH( (smooth_height/s2)**n -  ( (sum(dz_scl(1:(i-1)))+dz_scl(i)/2) /s2)**n ) / SINH((smooth_height/s2)**n)   ! small terrain features
  +                                + h1_u  *  SINH( (smooth_height/s1)**n -  ( (sum(dz_scl(1:(i-1)))+dz_scl(i)/2) /s1)**n ) / SINH((smooth_height/s1)**n)  &! large-scale terrain
  +                                + h2_u  *  SINH( (smooth_height/s2)**n -  ( (sum(dz_scl(1:(i-1)))+dz_scl(i)/2) /s2)**n ) / SINH((smooth_height/s2)**n)   ! small terrain features
                       z_v(:,i,:)   = (sum(dz_scl(1:(i-1))) + dz_scl(i)/2)   &
  -                                  ! + terrain_v  *  SINH( (H/s)**n - ( (sum(dz_scl(1:(i-1)))+dz_scl(i)/2) /s)**n ) / SINH((H/s)**n)
  -                                   + h1_v  *  SINH( (smooth_height/s1)**n -  ( (sum(dz_scl(1:(i-1)))+dz_scl(i)/2) /s1)**n ) / SINH((smooth_height/s1)**n)  &! large-scale terrain
  -                                   + h2_v  *  SINH( (smooth_height/s2)**n -  ( (sum(dz_scl(1:(i-1)))+dz_scl(i)/2) /s2)**n ) / SINH((smooth_height/s2)**n)   ! small terrain features
  +                                + h1_v  *  SINH( (smooth_height/s1)**n -  ( (sum(dz_scl(1:(i-1)))+dz_scl(i)/2) /s1)**n ) / SINH((smooth_height/s1)**n)  &! large-scale terrain
  +                                + h2_v  *  SINH( (smooth_height/s2)**n -  ( (sum(dz_scl(1:(i-1)))+dz_scl(i)/2) /s2)**n ) / SINH((smooth_height/s2)**n)   ! small terrain features
   
                       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 )
  
positive.code: 
      end subroutine initialize_core_variables
  
      subroutine setup_dzdxy(this,options)
          implicit none
          class(domain_t), intent(inout)  :: this
  
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
  
  