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:               endif
                dzfdx(:,i,:) = (zf(ims+1:ime,i,:) - zf(ims:ime-1,i,:)) / this%dx  
                dzfdy(:,i,:) = (zf(:,i,jms+1:jme) - zf(:,i,jms:jme-1)) / this%dx
                print*,"Level ",i,"  zf(150,i,:): ",zf(150,i,:)
              enddo
  
              ! Then finally:
  
negative.code:                                     this% kms : this% kme, &
                                      this% jms : this% jme) )
  
              print*,"maxval(forcing_terrain): ", maxval(forcing_terrain)     
              do i = this%grid%kms, this%grid%kme
                
                if (i<=max_level) then
  