anchor.diff:                  zf_interface(:,i+1,:)  = sum(dz_scl(1:i))   &
                                          + forcing_terrain  *  SINH( (H/s)**n - (sum(dz_scl(1:i))/s)**n ) / SINH((H/s)**n) 
   
  +                print*,"zf_inteface (150,150",i,"): ",zf_interface(:,i,:)
   
                 
                   if (i/=this%grid%kme)  dzf_interface(:,i,:)  =  zf_interface(:,i+1,:) - zf_interface(:,i,:) 
  
positive.code:             call io_write("delta_dzdx_sc.nc", "delta_dzdx", delta_dzdx(:,:,:) )
              call io_write("delta_dzdx_lc.nc", "delta_dzdx_lc", delta_dzdx_lc(:,:,:) )
              print*, "sinh term ", sin_h(:) 
              call io_write("dzdx.nc", "dzdx", dzdx(:,:) )
              call io_write("dzfdx.nc", "dzfdx", dzfdx(:,:) )
              call io_write("zf.nc", "zf", zf(:,:) )
  
      
          !_________ 2. Calculate the ratio bewteen z levels from hi-res and forcing data for wind acceleration  _________                                
  
negative.code:             ! h = terrain(:,:) 
  
              ! Scale dz with smooth_height/sum(dz(1:max_level)) before calculating sleve levels. 
              
              dz_scl(:)   =   dz(1:nz)  *  H / sum(dz(1:max_level))
  
              ! dz_scl   =   dz(:)  *  smooth_height / sum(dz(1:max_level))
              ! H        =  sum(dz_scl(1:max_level))  ! should also lead to smooth_height, but more error proof?
  