anchor.diff:              allocate(dzfdy(this% ims : this% ime, &
                             this% kms : this% kme, &
                             this% jms+1 : this% jme) )
  -
  -        ! allocate(this%dzdx(this% ims+1 : this% ime, &
  -        !                    this% kms : this% kme, &
  -        !                    this% jms : this% jme) )
  -
  -        ! allocate(this%dzdy(this% ims : this% ime, &
  -        !                    this% kms : this% kme, &
  -        !                    this% jms+1 : this% jme) )
  -
   
               allocate(delta_terrain(this% ims : this% ime, &
                                       this% jms : this% jme) )
  
positive.code: 
  
              ! Then finally:
              delta_dzdx(:,i,:) = dzdx(:,i,:)  -  dzfdx(:,i,:)  ! use this for w_real calculation.
  
              ! the same??? without all the above??? could be function in if(use_delta_terrain_for_speedup==True)
              delta_terrain = (terrain - forcing_terrain)
  
negative.code:             endif
  
  
              i=kms 
              
              ! - - - - -   Mass grid calculations for lowest level (i=kms)  - - - - -
              z_interface(:,i,:)   = terrain
  
              z_interface(:,i+1,:)  = dz_scl(i)   &
                                       + terrain  *  SINH( (H/s)**n - (dz_scl(i)/s)**n ) / SINH((H/s)**n)  ! same for higher k
                          
              dz_interface(:,i,:)  =  z_interface(:,i+1,:) - z_interface(:,i,:)  ! same for higher k
              
              z_level_ratio(:,i,:) = dz_interface(:,i,:)  / dz_scl(i)  ! same for higher k
              dz_mass(:,i,:)       = dz_interface(:,i,:) / 2           ! Diff for k=1
              z(:,i,:)             = terrain + dz_mass(:,i,:)          ! Diff for k=1   
  
              dzdx(:,i,:) = (z(ims+1:ime,i,:) - z(ims:ime-1,i,:)) / this%dx  ! same for higher k
              dzdy(:,i,:) = (z(:,i,jms+1:jme) - z(:,i,jms:jme-1)) / this%dx
  
  
              ! - - - - -   u/v grid calculations for lowest level (i=kms)  - - - - -
              ! for the u and v grids, z(1) was already initialized with terrain.
              ! but the first level needs to be offset, and the rest of the levels need to be created
              ! BK: So if z_u is already offset in the u dir, but not in the z dir, we can say that
              !     z_u(:,1,:) is the terrain on the u grid, and it needs to be offset in the z-dir 
              !     to reach mass levels (so by dz[i]/2)
  
              terrain_u =  z_u(:,i,:)  ! save for later on. 
              terrain_v =  z_v(:,i,:)  ! save for later on 
   
              ! 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)                                        
              z_v(:,i,:)  = dz_scl(i)/2   & 
                            + terrain_v  *  SINH( (H/s)**n - ( dz_scl(i)/2 /s)**n ) / SINH((H/s)**n)                                        
  
  
              zr_u(:,i,:)  =  (z_u(:,i,:) - terrain_u) / ( dz_scl(i)/2 )
              zr_v(:,i,:)  =  (z_v(:,i,:) - terrain_v) / (dz_scl(i)/2 )
  
              ! - - - - -  higher k levels  - - - - - 
              do i = this%grid%kms+1, this%grid%kme
  
                  if (i<=max_level) then  
  
                      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,:)  
                      
                      else  
  
                        z_interface(:,i+1,:)  =  sum(dz_scl(1:i))    &
                                       + terrain * SINH( (H/s)**n - (sum(dz_scl(1:i))/s)**n ) / SINH((H/s)**n)     
                      
                        dz_interface(:,i,:)  = z_interface(:,i+1,:) - z_interface(:,i,:)  
  
                      endif
                     
                      z_level_ratio(:,i,:) = dz_interface(:,i,:) / dz_scl(i)
  
                      ! if ( ANY(z_level_ratio(:,i,:)<0) ) then   ! Eror catching. Probably good to engage.
                      !   if (this_image()==1) then
                      !     write(*,*) "Error: z_level_ratio below zero (for level  ",i,")"
                      !     error stop
                      !     print*, dz_interface(:,i,:)
                      !     print*,""
                      !   endif
                      ! else if ( ANY(z_level_ratio(:,i,:)<=0.01) ) then
                      !   if (this_image()==1)  write(*,*) "WARNING: z_level_ratio 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)  
                      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)  
  
                      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 )
                  
                  else ! above the flat_z_height
                      
                      z_level_ratio(:,i,:) = 1
                      zr_u(:,i,:) = 1
                      zr_v(:,i,:) = 1
  
                      dz_interface(:,i,:) =   dz(i)
                      if (i/=this%grid%kme)   z_interface(:,i+1,:) = z_interface(:,i,:) + dz(i)
  
                      z_u(:,i,:)  = z_u(:,i-1,:)  + ((dz(i)/2 * zr_u(:,i,:) + dz(i-1)/2 * zr_u(:,i-1,:))) ! zr_u only relevant for first i above max level, aferwards both zr_u(i) AND zr_u(i-1) are 1
                      z_v(:,i,:)  = z_v(:,i-1,:)  + ((dz(i)/2 * zr_v(:,i,:) + dz(i-1)/2 * zr_v(:,i-1,:)))
  
                  endif
                  
                  dz_mass(:,i,:)   =  dz_interface(:,i-1,:) / 2  +  dz_interface(:,i,:) / 2
                  z(:,i,:)         =  z(:,i-1,:)           + dz_mass(:,i,:)
                  dzdx(:,i,:)      =  (z(ims+1:ime,i,:) - z(ims:ime-1,i,:)) / this%dx
                  dzdy(:,i,:)      =  (z(:,i,jms+1:jme) - z(:,i,jms:jme-1)) / this%dx
  
              enddo  ! ____ end SLEVE simple Implementation  _______
              
  
              ! if (options%parameters%debug) then 
              !   ! Write the z_level_ratio for analysis:
              !   write(file_idn, *) options%parameters%sleve_n
              !   write(file_idS, *) options%parameters%sleve_decay_factor
              !   filename = 'zlr_n_' // trim(adjustl(file_idn)) // '_S_'//trim(adjustl(file_idS)) //'.nc'
              !   call io_write(filename, "z_level_ratio", z_level_ratio(:,:,:) )
              !   ! call io_write("zlr.nc", "z_level_ratio", z_level_ratio(:,:,:) )
              !   ! call io_write("dzdx.nc", "dzdx", dzdx(:,:,:) )
              ! endif
  
  
            else  !. i.e. no sleve coordinates
              i = this%grid%kms
  
              max_level = nz
  
              if (options%parameters%space_varying_dz) then
                  max_level = find_flat_model_level(options, nz, dz)
  
                  smooth_height = sum(global_terrain) / size(global_terrain) + sum(dz(1:max_level))
  
                  z_level_ratio(:,i,:) = (smooth_height - terrain) / sum(dz(1:max_level))
                  global_z_level_ratio(:,i,:) = (smooth_height - global_terrain) / sum(dz(1:max_level))
  
                  zr_u(:,i,:) = (smooth_height - z_u(:,i,:)) / sum(dz(1:max_level))
                  zr_v(:,i,:) = (smooth_height - z_v(:,i,:)) / sum(dz(1:max_level))
              else
                  z_level_ratio = 1
                  global_z_level_ratio = 1
                  zr_u = 1
                  zr_v = 1
              endif
  
              dz_mass(:,i,:)      = dz(i) / 2 * z_level_ratio(:,i,:)
              dz_interface(:,i,:) = dz(i) * z_level_ratio(:,i,:)
              z(:,i,:)            = terrain + dz_mass(:,i,:)
              z_interface(:,i,:)  = terrain
  
              global_dz_interface(:,i,:) = dz(i) * global_z_level_ratio(:,i,:)
              global_z_interface(:,i,:)  = global_terrain
  
  
              ! for the u and v grids, z(1) was already initialized with terrain.
              ! but the first level needs to be offset, and the rest of the levels need to be created
              z_u(:,i,:)          = z_u(:,i,:) + dz(i) / 2 * zr_u(:,i,:)
              z_v(:,i,:)          = z_v(:,i,:) + dz(i) / 2 * zr_v(:,i,:)
  
              dzdx(:,i,:) = (z(ims+1:ime,i,:) - z(ims:ime-1,i,:)) / this%dx
              dzdy(:,i,:) = (z(:,i,jms+1:jme) - z(:,i,jms:jme-1)) / this%dx
  
              do i = this%grid%kms+1, this%grid%kme
                  if (i<=max_level) then
                      z_level_ratio(:,i,:) = z_level_ratio(:,i-1,:)
                      zr_u(:,i,:) = zr_u(:,i-1,:)
                      zr_v(:,i,:) = zr_v(:,i-1,:)
  
                      global_z_level_ratio(:,i,:) = global_z_level_ratio(:,i-1,:)
  
                  else
                      z_level_ratio(:,i,:) = 1
                      zr_u(:,i,:) = 1
                      zr_v(:,i,:) = 1
  
                      global_z_level_ratio(:,i,:) = 1
  
                  endif
  
                  dz_mass(:,i,:)     = (dz(i)/2 * z_level_ratio(:,i,:) + dz(i-1)/2 * z_level_ratio(:,i-1,:))
                  dz_interface(:,i,:)= dz(i) * z_level_ratio(:,i,:)
                  z(:,i,:)           = z(:,i-1,:)           + dz_mass(:,i,:)
                  z_interface(:,i,:) = z_interface(:,i-1,:) + dz_interface(:,i-1,:)
  
                  global_dz_interface(:,i,:) = dz(i) * global_z_level_ratio(:,i,:)
                  global_z_interface(:,i,:)  = global_z_interface(:,i-1,:) + global_dz_interface(:,i-1,:)
  
                  z_u(:,i,:)         = z_u(:,i-1,:)         + ((dz(i)/2 * zr_u(:,i,:) + dz(i-1)/2 * zr_u(:,i-1,:)))
                  z_v(:,i,:)         = z_v(:,i-1,:)         + ((dz(i)/2 * zr_v(:,i,:) + dz(i-1)/2 * zr_v(:,i-1,:)))
  
                  dzdx(:,i,:) = (z(ims+1:ime,i,:) - z(ims:ime-1,i,:)) / this%dx
                  dzdy(:,i,:) = (z(:,i,jms+1:jme) - z(:,i,jms:jme-1)) / this%dx
              enddo
   
          endif
  
              i = this%grid%kme + 1
              global_z_interface(:,i,:) = global_z_interface(:,i-1,:) + global_dz_interface(:,i-1,:)
  
              ! 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,:)
              ! dz_mass(:,i,:)     = dz(i-1)/2 * z_level_ratio(:,i-1,:)
  
          use_delta_terrain=.true.    
          if(use_delta_terrain) then
          
             ! #----------------------- calc z levels from forcing terrain -------------------
              print* , " hgt variable forcing_terrain: ", shape(this%forcing_terrain%data_2d)
              call io_write("forcing_terrain.nc", "forcing_terrain", forcing_terrain(:,:) ) ! check in plot
  
              ! - - - - - - - - - - - - - -
              allocate( zf_interface(this% ims : this% ime, &
                                      this% kms : this% kme+1, &
                                      this% jms : this% jme) )
  
              allocate(dzf_interface(this% ims : this% ime, &
  