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:               if (i<=max_level) then
                
                  if (i==this%grid%kms)    zf_interface(:,i,:)   =  forcing_terrain
                  if (i==this%grid%kme)    dzf_interface(:,i,:)  =  smooth_height - zf_interface(:,i,:)  
                 
                  zf_interface(:,i+1,:)  = dz_scl(i)   &
                                         + forcing_terrain  *  SINH( (H/s)**n - (dz_scl(i)/s)**n ) / SINH((H/s)**n)  ! same for higher k
  
negative.code:             ! dz_mass(:,i,:)     = dz(i-1)/2 * z_level_ratio(:,i-1,:)
  
      
          end associate
  
          temp =  this%zr_u
  