anchor.diff:              ! get the associated forcing data
               input_data = forcing%variables%get_var(var_to_interpolate%forcing_var)
   
  -            if (var_to_interpolate=="forcing_terrain") print*,"interpolating forcing_terrain"
  -
  +            
               ! interpolate
               if (var_to_interpolate%two_d) then
                   if (update_only) then
                       call geo_interp2d(var_to_interpolate%dqdt_2d, input_data%data_2d, forcing%geo%geolut)
                   else
                       call geo_interp2d(var_to_interpolate%data_2d, input_data%data_2d, forcing%geo%geolut)
  -                    if (var_to_interpolate=="HSURF") print*,"interpolated forcing_terrain ", shape(input_data%data_2d)
                   endif
   
               else
  
positive.code:                   terrain_v             => this%terrain_v%data_2d,              &
                    forcing_terrain       => this%forcing_terrain%data_2d,        &
                    n                     =>  options%parameters%sleve_n,         & 
                    dzdx                  => this%dzdx,                           &
                    dzdy                  => this%dzdy,                           &
                    dz_scl                => this%dz_scl,                          &
  
negative.code:             ! dz_mass(:,i,:)     = dz(i-1)/2 * z_level_ratio(:,i-1,:)
  
      
          end associate
  
          temp =  this%zr_u
  