anchor.diff:                  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
  +
  +
  +            ! # # # # # #   BK debug wind field:    # # # # # # # # #
  +            if (options%parameters%debug ) then
  +                print*, ""
  +                print*, "z_level_ratio(150,i,150) : ", z_level_ratio(150,:,150)
  +                print*, "dzdx at  i,j,k = 150,150,k: " !, k, " : "
  +                print*, dzdx(150,:,150)
  +                print*, "dzdy at i,j,k = 150,150,k: "!, k, " : "
  +                print*, dzdy(150,:,150)
  +
  +            endif  
  +
   
   
           end associate
  
positive.code: 
          if (this_image()==1) print *,"  Initializing variables"
  
          if (0<opt%vars_to_allocate( kVARS%u) )                          call setup(this%u,                        this%u_grid,   forcing_var=opt%parameters%uvar,       list=this%variables_to_force, force_boundaries=.False.)
          if (0<opt%vars_to_allocate( kVARS%u) )                          call setup(this%u_mass,                   this%grid)
          if (0<opt%vars_to_allocate( kVARS%v) )                          call setup(this%v,                        this%v_grid,   forcing_var=opt%parameters%vvar,       list=this%variables_to_force, force_boundaries=.False.)
          if (0<opt%vars_to_allocate( kVARS%v) )                          call setup(this%v_mass,                   this%grid)
          if (0<opt%vars_to_allocate( kVARS%w) )                          call setup(this%w,                        this%grid )
          if (0<opt%vars_to_allocate( kVARS%w) )                          call setup(this%w_real,                   this%grid )
          if (0<opt%vars_to_allocate( kVARS%water_vapor) )                call setup(this%water_vapor,              this%grid,     forcing_var=opt%parameters%qvvar,      list=this%variables_to_force, force_boundaries=.True.)
          if (0<opt%vars_to_allocate( kVARS%potential_temperature) )      call setup(this%potential_temperature,    this%grid,     forcing_var=opt%parameters%tvar,       list=this%variables_to_force, force_boundaries=.True.)
          if (0<opt%vars_to_allocate( kVARS%cloud_water) )                call setup(this%cloud_water_mass,         this%grid,     forcing_var=opt%parameters%qcvar,      list=this%variables_to_force, force_boundaries=.True.)
  
negative.code: 
                  smooth_height = sum(this%global_terrain) / size(this%global_terrain) + sum(dz(1:max_level))
  
                  z_level_ratio(:,i,:) = (smooth_height - terrain) / sum(dz(1:max_level)) !BK: Ratio between spatially varying and non-varying dz(at flatdz height) 
  
                  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))
  