anchor.diff:                                      this% kms : this% kme, &
                                       this% jms : this% jme) )
   
  -            print*,"maxval(forcing_terrain): " maxval(forcing_terrain)     
  +            print*,"maxval(forcing_terrain): ", maxval(forcing_terrain)     
               do i = this%grid%kms, this%grid%kme
                 
                 if (i<=max_level) then
  
positive.code:             ! h = terrain(:,:) 
  
              ! Scale dz with smooth_height/sum(dz(1:max_level)) before calculating sleve levels. 
              print *, "dz shape: ", shape(dz)
              dz_scl(:)   =   dz(:)  *  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?
  
negative.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  _________                                
  