anchor.diff:              ! #------------------------------------------
   
   
  -            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,:)
  +        endif    
  +
   
           end associate
   
  
positive.code:             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) )
  
negative.code: 
        this%model_time = forcing%current_time
  
        ! call terrain_delta(this, forcing, options)  ! BK TEST
      end subroutine
  
  
  